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I.  Summary 


The  research  performed  under  the  contract,  during  the  period  14  July  1985  through 
13  January  1986,  can  be  divided  into  two  main  topics;  Ahe  effects  of  lateral  structure,  in 
particular,  complicated  sedimentary  basins  on  the  propagation  and  amplication  of 
seismic  strong  motion;  andjthe  effects  of  ocean  continent  transition  zones  on  Lt  waves.  1) 

In  section  II,  2-D  finite  diffferpce  techniques  are  used  to  model  the  seismic  strong 
motions  observed  for  the  1971  San  Fernando  earthquake.  The  finite  difference  seismo¬ 
grams  show  strong  effects  due  to  lateral  variation  in  sediment  thickness  in  the  San  Fer¬ 
nando  valley  and  the  Los  Angeles  basin.  Using  basin  structure  derived  mostly  from  well 
logs  and  teleseismically  determined  source  parameters,  two  dimensional  SH  and  P-SV 
finite  difference  calculations  reproduce  the  amplitude  and  duration  of  the  strong  motion 
velocities  recorded  across  the  Los  Angeles  and  San  Fernando  basins  for  the  period  range 
2  to  10  seconds.  The  edges  of  basins  nearest  the  seismic  source  show  ground  motion 
amplification  up  to  a  factor  of  three  over  the  case  without  the  basin,  and  tend  to  convert 
direct  shear  waves  into  Love  and  Rayleigh  waves  that  travel  within  the  basin.  )  > 

In  section  III,  the  methods  for  Representation  Theorem^RT)"coupling  of  finite  ele¬ 
ment  {FE)-or  finite  difference  calculations  and  propagator  matrix  method  calculations  to 
produce  a  hybrid  method  for  propagation  of  SH  mode  sum  seismograms  across  paths 
that  contain  regions  of  non  plane-layered  structure  are  explained  and  developed  AThe 
coupling  methods  explained  in  detail  use  a  2-D  Cartesian  FE  formulation.  Extensive 
tests  illustrating  the  validity  and  accuracy  of  the  implementation  of  these  coupling 
methods  are  discussed.  These  hybrid  techniques  are  developed  to  study  the  propagation 
of  surface  waves  across  regional  transition  zones  or  other  heterogeneities  that  exist  in 
part  of  a  longer,  mostly  plane-layered,  path.  The  effects  of  a  thinning  or  thickening  of 
the  crustal  layer  on  the  propagation  of  Lt  mode  sum  seismograms  will  be  examined  in 
later  reports. 
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ABSTRACT 


Finite-difference  seismograms  calculated  for  the  1971  San  Fernando  earth¬ 
quake  show  strong  effects  due  to  lateral  variation  in  sediment  thickness  in  the 
San  Fernando  valley  and  the  Los  Angeles  basin.  Using  basin  structure  derived 
mostly  from  well  logs  and  teleseismically-determined  source  parameters,  two- 
dimensional  SH  and  P-SV  finite  difference  calculations  reproduce  the  amplitude 
and  duration  of  the  strong  motion  velocities  recorded  across  the  Los  Angeles  and 
San  Fernando  basins  for  the  period  range  2  to  10  seconds.  The  edges  of  basins 
nearest  the  seismic  source  show  ground  motion  amplification  up  to  a  factor  of 
three  over  the  case  without  the  basin,  and  tend  to  convert  direct  shear  waves 
into  Love  and  Rayleigh  waves  that  travel  within  the  basins. 


INTRODUCTION 


The  San  Fernando  earthquake  of  February  9,  1971  produced  one  of  the  most 
complete  sets  of  strong-ground  records  from  a  large  thrust  earthquake  to  date, 
and  consequently  has  generated  a  large  body  of  seismological  literature.  Hanks 
(1975)  notes  that  seismic  moment,  source  dimension,  radiation  pattern,  rupture 
propagation,  the  development  and  dispersion  of  local  surface  waves,  and  azimu¬ 
thal  variations  in  the  gross  geological  structure  all  appear  to  have  first-order 
significance  in  fashioning  the  amplitude  and  frequency  content  of  the  strong- 
motion  waveforms.  Heaton  (1982),  among  others,  has  modeled  some  of  these 
features,  but  adds  that  many  features  of  the  observed  motions  remain  unex¬ 
plained,  and  considerable  uncertainty  still  exists  regarding  the  faulting  history  of 
the  San  Fernando  earthquake. 

One  drawback  of  existing  studies  is  the  inability  to  properly  simulate  the 
propagational  effects  associated  with  the  complicated  sedimentary  basins  in  and 
around  Los  Angeles.  Finite  element  schemes  have  been  used  to  explain  some  of 
the  effect  of  basins  on  surface  waves  body  waves,  and  the  static  deformation 
(Drake  and  Mai,  1972,  McCowan  et  al.,  1977).  However,  these  efforts  did  not 
address  the  demanding  task  of  modeling  seismic  profiles  along  paths  crossing  the 
various  basin  and  ridge  structures,  as  discussed  qualitatively  by  Liu  and  Heaton 
(1984).  The  modeling  of  one  such  three-component  profile  is  the  subject  of  this 
study. 

STRONG  MOTION  RECORDS 

The  filled  triangles  in  Figure  1  locate  three  component  accelerometers  that 
recorded  the  strong  motions  of  the  San  Fernando  earthquake  along  a  north-south 
profile.  Figure  2  shows  the  vertical,  radial,  and  transverse  velocities  recorded  at 
these  stations.  The  accelerations  were  integrated  to  velocities  by  EERL  (1974). 


We  model  velocities  two  reasons;  first,  it  is  currently  economical  to  model 
only  greater  than  1  sec  period  seismic  energy  to  ranges  of  50-60  km  with  our 
finite-difference  scheme,  and  velocity  but  not  acceleration  records  have  these 
periods,  and  second,  we  have  less  confidence  in  our  knowledge  of  the  geologic 
structures  with  wavelengths  shorter  than  a  few  kilometers,  which  primarily 
influence  shorter-period  seismic  waves.  These  records  comprise  profile  I  in  Liu 
and  Heaton  (1984).  The  absolute  timing  of  these  records  is  not  known  so  they 
are  lined  up  relative  to  an  early,  high-frequency  arrival  on  the  vertical  component 
that  is  probably  a  direct  compressional  wave. 

Before  discussing  these  motions  in  detail,  it  is  useful  to  review  the  geological 
structure  along  this  profile.  A  schematic  cross-section  is  given  on  the  bottom  of 
Figure  1.  Figure  3  shows  the  geologic  layering  in  this  cross-section  in  more 
detail,  with  an  inferred  shear  wave  velocity  structure  below.  The  stratigraphy  is 
taken  from  Duke  et  al.  (1971),  who  reviewed  the  well  logs  and  the  geological 
cross-sections  in  the  literature  for  this  area.  Many  well  logs  have  been  recorded 
around  the  Los  Angeles  basins  because  of  oil  exploration  and  these  logs  yield  esti¬ 
mates  of  density  and  velocity.  Duke  et  al.  (1971)  also  conducted  numerous 
small-scale  refraction  surveys  to  find  the  near-surface  compressional-wave  velocity 
profile.  With  knowledge  of  the  composition  of  the  rocks  at  depth,  the  shear-wave 
velocity  can  also  be  reliably  estimated.  Near  surface  values  of  the  attenuation 
operator  Q  have  also  been  provided  by  a  variety  of  schemes.  Duke  et  al.  (1971) 
reports  the  details  of  this  cross-section  at  great  length. 

Although  the  structure  within  the  San  Fernando  and  Los  Angeles  basins 
along  this  profile  can  be  reliably  estimated  from  the  above  study,  the  structure 
below  the  basins  is  less  well-known  since,  unfortunately,  the  well  logs  stop  as 
soon  as  they  encounter  the  basement  rock.  Below  the  basins,  we  use  the  struc¬ 
ture  given  in  Kanamori  and  Hadley  (1975)  for  compressional  waves,  and  assume 


that  Vp  /  Vt  is  s/l.  Table- 1  gives  the  compressional  wave  velocity  and  density 
associated  with  each  shear  wave  velocity  shown  in  Figure  3. 

The  motions  shown  in  Figure  2  correspond  with  the  geologic  setting  in  which 
they  were  recorded.  Within  the  San  Fernando  basin,  a  train  of  surface  waves 
develops  with  an  apparent  velocity  of  1  km/sec.  The  wavetrain  lengthens  as  it 
propagates  across  the  basin.  On  the  Santa  Monica  mountains,  the  surface  waves 
disappear  and  the  amplitude  falls  by  a  factor  of  two.  In  the  Los  Angeles  basin, 
the  surface  wave  is  again  present  with  a  slow  apparent  velocity,  and  the  ampli¬ 
tudes  increase  by  roughly  a  factor  of  two  relative  to  the  ridge.  Near  the  ocean, 
in  the  Palos  Verdes  hills,  the  surface  wave  is  still  present,  but  the  amplitude  has 
dropped  more  rapidly  than  geometric  spreading  would  predict.  Geometric 
spreading  would  be  1  f\fr  for  surface  waves  and  1/r  for  body  waves  in  a 
wholespace.  Spreading  for  a  layered  halfspace  can  fall  off  by  higher  powers  of  r 
for  body  waves,  but  the  data  suggest  surface  waves  are  present  and  show  a 
decrease  by  a  factor  of  3  from  Hi  18  and  S267  to  N191,  where  the  range  increases 
from  50  to  70  km,  which  would  require  a  large  exponent  on  r.  It  is  these  pat¬ 
terns,  some  of  which  are  discussed  by  Liu  and  Heaton  (1984),  that  we  will 
attempt  to  understand  by  forward  modeling  through  the  FD  technique  with 
laterally  varying  structures. 

Before  forward  modeling  with  a  technique  that  assumes  two-dimensional 
symmetry,  we  will  verify  that  the  seismic  energy  in  the  records  is  not  significantly 
laterally  reflected.  Vidale  (1986)  uses  complex  polarization  analysis  that  is  an 
extension  of  the  method  of  Montalbetti  and  Kanasewich  (1970)  to  examine  the 
three-component  records  for  this  profile  of  records  from  the  San  Fernando  earth¬ 
quake.  Because  the  method  works  in  the  complex  domain,  elliptically  as  well  as 
linearly  polarized  data  may  be  interpreted.  This  analysis  finds  the  direction  of 
propagation  for  the  Love  and  Rayleigh  waves.  The  direction  may  be  obtained 


because  the  Love  wave  is  linearly  polarized  transverse  to  the  direction  of  propa¬ 
gation,  and  the  Rayleigh  wave  is  elliptically  polarized  in  the  plane  that  contains 
the  vertical  and  the  propagation  directions. 

The  analysis  in  Vidale  (1986)  shows  that  the  energy  in  this  particular  profile 
is  traveling  within  15°  of  radially  outward  from  the  source.  This  suggests  that 
despite  the  three-dimensional  nature  of  the  basins,  the  geometry  may  be  approxi¬ 
mated  by  a  two-dimensional  model  with  useful  results.  We  should  note,  however, 
that  although  the  energy  in  this  profile  travels  radially  out  from  the  source  to  the 
receiver,  amplitude  attenuation  due  to  geometrical  spreading  depends  on  the  cur¬ 
vature  of  the  wavefront,  so  the  amplitude  may  vary  somewhat  due  to  three- 
dimensional  effects  even  if  the  energy  path  is  not  laterally  deflected. 

To  match  the  frequency  limitations  of  the  FD  algorithm,  we  low-pass  filter 
the  data.  The  FD  method  can  only  propagate  energy  with  six  or  more  grid 
points  per  wavelength,  so  to  properly  treat  high  frequencies,  more  grid  points  and 
therefore  more  calculation  is  required.  Figure  2b  shows  the  filtered  vertical, 
radial,  and  transverse  velocities  that  will  be  addressed  in  this  modeling  study. 

REVIEW  OF  SOURCE  PARAMETERS 

The  San  Fernando  earthquake  has  also  been  the  subject  of  studies  that  con¬ 
centrated  on  other  data  sets,  namely  the  teleseismic  body  waves,  the  long-period 
surface  waves,  the  location  of  preshocks  and  aftershocks,  and  the  static  displace¬ 
ments.  These  studies  will  help  us  estimate  the  source  parameters.  First,  we  will 
briefly  review  these  studies  to  help  clarify  some  of  the  uncertainties. 

The  seismicity  studies  (Allen  et  al.,  1973,  Whitcomb  et  al.,  1973)  reveal  a 
relatively  diffuse  zone  of  aftershocks  with  a  combination  of  thrust  and  left-lateral 
strike-siip  mechanisms.  By  locating  the  mainshock  relative  to  well-located  aft¬ 
ershocks  recorded  on  a  temporarily  deployed  array,  Whitcomb  et  al.  (1973)  place 
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the  hypocenter  :.t  a  depth  of  8  km  at  34°  24.7'  N  and  118°  24.C/  W,  with  the 
mechanism  listed  in  Table  2.  The  errors  in  hypocentral  depth  are  conservatively 
estimated  at  ±  8  km  vertically  and  ±  4  km  horizontally.  The  aftershocks  sug¬ 
gest  that  the  plane  of  faulting  dipped  35°  to  the  north  down  to  8  km  depth, 
below  which  the  fault  plane  dipped  50° .  The  mechanism  derived  from  first 
motions  by  Whitcomb  et  al.  (1973),  listed  in  Table  2,  is  similar  to  the  teleseismic 
results. 

Examination  of  the  static  displacements  due  to  the  earthquake  (Alewine, 
1974)  reveals  distributed  slip  from  near  the  surface  to  a  depth  of  about  14  km 
along  the  north  dipping  fault  plane.  The  portion  of  the  fault  from  the  surface 
down  to  5  km  depth  underwent  about  5  m  of  slip,  and  the  segment  from  10  to  14 
km  depth  shows  2-5  m  slip,  and  there  may  be  an  area  of  less  slip  from  5  to  10  km 
depth.  The  greater  the  depth,  however,  the  worse  the  resolution  of  static 
analysis.  The  static  moment  is  estimated  to  be  between  1.0  and  2.2  X  1026 
dyne-cm. 

The  16-60  sec  fundamental  Rayleigh  waves  generated  by  the  San  Fernando 
earthquake  are  analyzed  by  Alewine  (1974),  and  given  the  strike  and  dip  sug¬ 
gested  by  Whitcomb  et  al.  (1973),  the  moment  is  found  to  be  1.7  X  1026  dyne-cm. 
Half  the  moment  release  is  found  to  be  at  depths  of  about  3-8  km  and  the  other 
half  of  the  moment  release  is  below  about  10  km. 

The  teleseismic  body  wave  studies  (Langston,  1978,  and  Heaton,  1982)  use 
records  from  long-period  WWSSN  stations  to  find  the  mechanism  and  5-30 
second  period  faulting  history  in  time  and  space. 

Both  Langston  (1978)  and  Heaton  (1982)  find  a  double  source.  In  the  vari¬ 
ous  models,  there  is  a  source  at  10-15  km  depth  with  a  moment  of  0.5  X  1026 
dyne-cm  and  there  is  a  shallower  source  with  a  shallower  dip  angle  that  also  has 
a  moment  of  0.5  X  1026.  The  two  sources  are  found  to  have  a  slightly  different 


strikes  and  the  shallower  source  dips  less,  as  shown  in  Table  2.  When  Heaton 
(1082)  attempts  to  model  the  near-in  records  and  static  data  as  well  as  the  telese- 
ismic  data,  his  estimate  of  the  total  moment  rises  to  1.7  X  1026  dyne-cm.  Figure 
4  summarizes  the  results  from  the  teleseismic  modeling  of  Langston  (1978)  and 
Heaton  (1982). 

The  short-period  WWSSN  records  are  also  examined  in  the  search  for  higher 
frequency  (0.5-3.0  second  period)  details  of  the  faulting  history  (Hanks,  1975,  and 
Langston,  1978).  The  work  of  both  Langston  (1978)  and  Hanks  (1975)  suggest 
that  the  first  pulse  of  short-period  energy  originated  about  12  km  below  the  sur¬ 
face.  If  the  pulse  of  energy  came  from  the  hypocenter,  which  was  on  the  fault 
plane  as  defined  by  the  aftershocks,  the  location  from  Whitcomb  et  al.  (1973) 
would  have  to  be  in  error  by  4  km  vertically  and  4  km  horizontally.  Some  of  the 
short-period  WWSSN  records,  however,  indicate  a  small  precursor  about  2 
seconds  before  the  initial  large  pulse,  so  perhaps  the  hypocenter  and  the  location 
of  the  short-period  energy  release  are  not  coincident. 

In  this  study,  we  find  that  a  point  source  at  10  km  depth  with  the  mechan¬ 
ism  of  Heaton’s  deep  source  (which  is  nearly  identical  to  Langston’s  deep  source) 
can  explain  most  of  the  data.  The  sensitivity  of  our  modeling  to  changes  in 
source  depth  and  mechanism  are  discussed  below. 

NUMERICAL  MODELS  OF  THE  STRONG  MOTIONS 

Our  methods  are  described  in  the  article  in  this  issue  by  Helmberger  and 
Vidale  (1987).  More  details  of  these  two-dimensional  finite-difference  (FD)  SH 
and  P-SV  algorithms  that  are  used  to  calculate  the  Green's  functions  are 
described  in  Vidale  et  al.  (1985),  Vidale  and  Helmberger  (1987),  Clayton  and 
Vidale  (1987),  and  Vidale  and  Clayton  (1986).  The  finite-difference  method  is  of 
fourth-order  accuracy  as  defined  by  Alford  et  al.  (1974).  The  seismic  energy  from 


the  source  is  introduced  into  the  grid  by  the  method  of  Alterman  and  Karal 
(1968).  This  method,  as  we  use  it,  requires  the  source  to  be  placed  in  a  region  of 
the  grid  of  constant  density  and  velocities. 

Expressions  describing  line-source  excitation  functions  used  in  the  line-to- 
point-source  mapping  may  be  found  in  Helmberger  and  Vidale  (1987).  The  exact 
dislocation  solution  is  expressed  in  an  asymptotic  series  in  which  the  displace¬ 
ment  decouples  into  the  SH  and  P-SV  systems  where  the  vertical  and  azimuthal 
radiation  patterns  separate  explicitly.  Therefore  we  can  make  separate  FD  calcu¬ 
lations  to  model  the  SH  and  P-SV  radiation.  The  expressions  containing  the 
vertical  radiation  patterns  are  loaded  into  the  grid  source  box.  The  azimuthal 
radiation  patterns  are  factored  in  later  according  to  the  fault  orientation.  This 
FD  procedure  does  not  handle  the  near-field  terms  properly  at  small  horizontal 
distances,  but  numerical  checks  indicate  accurate  results  beyond  about  10  km  for 
a  source  at  a  depth  of  10  km. 

The  Green’s  functions  for  the  transverse  component  of  velocity  for  the  four 
source  locations  shown  in  Figure  3  are  displayed  in  Figures  5  and  6  where  j  —4 
corresponds  to  the  strike-slip  and  j  =5  to  the  dip-slip  components.  The  Green’s 
functions  correspond  to  the  Vj  in  Equations  1  and  14  of  Vidale  et  al.  (1985),  and 
are  equivalent  to  Equation  24  in  Helmberger  and  Vidale  (1987,  this  issue)  for  the 
case  where  the  source  is  impulsive  and  only  one  of  A  4  and  A  5  is  nonzero.  A  Q  of 
25  is  included  in  the  calculation. 

The  ranges  indicated  on  the  left  are  appropriate  for  the  deepest  source.  For 
shallower  sources,  we  chose  to  leave  the  receivers  at  the  same  point  in  the  struc¬ 
ture  with  the  same  range,  rather  than  shift  them  by  the  difference  in  hypocentral 
position,  so  the  ranges  are  approximate  for  the  shallower  sources.  Note  that  the 
shallower  sources  are  located  further  to  the  south.  The  change  in  polarity  for  the 
first  station  for  the  source  S  indicates  that  the  source  is  just  to  the  south  of  the 


first  receiver.  However,  we  are  only  interested  in  the  records  written  at  the  larger 
distances  where  the  near-field  terms  can  be  neglected.  A  seismic  attenuation  Q  of 
25  is  assumed;  this  value  falls  within  the  wide  range  of  Q’s  reported  by  Duke  et 
al.  (1971). 

Incorporation  of  spatially  varying  attenuation  into  finite-difference  calcula¬ 
tions  'is  possible  with  Pade  approximants  (Day  and  Minster,  1984,  Witte  and 
Richards,  1986)  but  requires  much  computation.  We  use  a  simpler  scheme  in  this 
paper  to  include  the  effect  of  a  constant  Q.  Our  scheme  is  derived  from  one  that 
has  been  used  in  generalized  ray  calculations  for  some  years  (Carpenter,  1967). 
The  quantity  t*  =  Jl/vQ  ds  or  t*  =  Jl/Q  dt  may  be  computed  for  a  given 
ray,  and  a  convolutional  attenuation  operator  can  be  computed  for  any  t  *  (see 
Futterman,  1962,  for  example).  If  Q  is  constant,  t  *  =  t  /Q ,  thus  t*  is  only  a 
function  of  time.  In  the  computation  of  teleseismic  arrivals,  t  *  remains  approxi¬ 
mately  constant  during  the  time  window  in  which  a  P-wave  or  even  several 
branches  of  a  P-wave  triplication  arrive. 

The  profile  in  this  paper  is  more  complicated,  however,  in  that  the  time  win¬ 
dow  is  not  short  compared  to  the  travel  time  from  the  source.  Since  Q  is 
included  with  a  convolution-like  operator,  however,  we  need  not  use  the  same 
operator  for  each  point  in  time.  Q  is  then  incorporated  by  convolving  at  each 
point  in  time  of  the  seismograms  the  Futterman  Q  operator  for  a  t*  ==  t/Q. 
Of  course,  the  convolution-like  operation  may  not  be  done  in  the  frequency 
domain  since  every  time  point  requires  a  different  Q  operator.  In  practice,  we 
only  calculate  a  new  Q  operator  for  each  50  time  points  since  the  Q  operator 
changes  gradually  with  time,  and  computing  an  operator  requires  an  FFT. 

A  moment  of  1026  dyne-cm  is  used  for  the  10  km  source  depth.  To  compare 
amplitudes  between  different  depths,  we  hold  AD0  rather  than  the  moment 
nAD0  constant  to  avoid  the  strong  tendency  to  generate  larger  seismic  motions 


in  softer 


The  Love  waves  dominate  the  motions  for  the  dip-slip  case  at  all  depths 
whereas  the  direct  S  arrivals  are  more  noticeable  in  the  strike-slip  case,  especially 
at  the  greater  depths.  These  features  may  be  understood  in  terms  of  vertical 
radiation  patterns.  The  strike-slip  source  radiates  most  of  the  energy  horizon¬ 
tally,  which  appears  as  direct  body  waves,  while  the  dip-slip  pattern  tends  to 
radiate  energy  vertically,  where  it  can  be  trapped  to  form  surface  waves. 

The  strike-slip  Green’s  function  is  the  most  important  for  the  transverse 
component  of  motion  for  the  San  Fernando  earthquake,  despite  the  fact  that  the 
mechanism  is  dominantly  thrust.  This  is  because  a  45°  dip-slip  event  produces 
the  same  radiation  pattern  of  SH  energy  as  a  strike-slip  event,  but  rotated  45°  in 
strike.  The  A  4  coefficient  (strike-slip  case)  is  about  10  times  larger  than  the  A  5 
coefficient  (dip-slip  case)  in  both  Langston's  (1078)  and  Heaton’s  (1082)  solution 
for  the  deep  source.  Heaton’s  shallow  source  is  also  dominated  by  the  strike-slip 
solution  but  Langston’s  shallow  source  is  about  one-third  composed  of  the  dip- 
slip  solution.  For  these  reasons,  we  have  conducted  most  of  our  sensitivity  stu¬ 
dies  with  the  strike-slip  case.  Synthetic  transverse  ground  motions  can  be  gen¬ 
erated  from  the  Green’s  functions  shown  in  Figures  5  and  6  by  performing  the 
operations  indicated  in  equation  1  of  Vidale  et  al.  (1985).  Transverse  velocity 
records  convolved  with  the  3  second  time  function  shown  in  Figure  8c  are  shown 
in  Figure  7b.  Seismograms  for  a  flat-layered  approximation  to  the  deep  basin 
structure  are  included  in  Figure  7a  for  comparison. 

These  synthetic  seismograms  for  a  flat  structure  are  sensitive  to  the  source 
depth  in  both  amplitude  and  waveform.  The  shallowest  source  excites  large,  slow 
Love  waves,  and  would  make  even  larger,  higher-frequency  waves  were  it  not  for 
the  strong  damping  due  to  the  Q  of  25.  Even  considering  the  inefficiency  of  the 
mechanism  of  the  shallow  source  in  exciting  radiation  along  this  profile,  0.1  X 


1026  dyne-cm  of  moment  in  the  shallow  source  would  suffice  to  generate  waves  as 
large  as  those  observed. 

The  shallow  source  creates  such  large  waves  for  two  reasons.  First,  the 
softer  material  near  the  surface  allows  larger  velocities  to  develop  for  a  given 
moment  than  the  stiffer  material  at  greater  depth.  Second,  the  shallower  source 
is  closer  to  the  edge  of  the  slow  San  Fernando  basin  than  the  deeper  sources,  and 
it  can  trap  a  larger  portion  of  its  energy  as  surface  waves.  The  direct  body-wave 
has  a  negligible  amplitude  compared  to  the  surface  wave  for  all  but  the  two 
closest  stations.  Source  I  generates  smaller  seismic  waves,  and  excites  surface 
waves  with  a  range  of  velocities.  The  direct  arrival  is  distinct  for  all  ranges,  and 
source  D  produces  predominantly  a  direct  diving  wave.  Surface  waves  do  not 
contribute  much  energy  to  the  records  for  the  intermediate  and  deep  sources.  It 
is  clear  that  for  this  flat-layered  geometry,  the  shallow  source  generates  much 
more  surface  motion  per  unit  moment  than  the  intermediate  and  deep  sources. 

The  moments  of  1-2  X  1026  dyne-cm  found  by  the  studies  in  Table  2,  placed 
at  the  intermediate  or  deep  depth  are  consistant  with  the  peak  velocities 
observed.  The  synthetic  waveform  fit  with  the  data,  however,  is  improved  con¬ 
siderably  by  the  empirical  source  time  function  described  below. 

The  waveforms  from  the  more  realistic  structure,  Figure  7b,  are  more  com¬ 
plicated  than  in  the  flat-layered  case.  Sources  at  all  depths  can  excite  noticeable 
surface  waves  in  both  basins.  This  result  arises  because  direct  energy  incident  on 
the  edges  of  the  basin  can  be  trapped  by  the  basins  much  more  efficiently  than 
the  direct  waves  can  tunnel  into  the  flat  layers  (Vidale  et  al.,  1985,  for  example). 
The  source  S  still  excites  larger  motions  on  the  surface,  again  largely  because  a 
source  in  the  slower  medium  generates  larger  amplitude  waves,  and  also  because 
the  shallow  source  is  closer  to  the  basin  edge,  allowing  a  larger  percentage  of  its 
energy  to  be  trapped.  A  source  at  location  B  produces  seismograms  very  similar 


to  those  of  a  s^jive  at  location  D,  except  for  a  smaller  amplitude,  which  is  due  to 
the  faster  velocity  at  the  deeper  depth. 

The  surface  waves  that  are  traveling  within  the  San  Fernando  basin  can,  to 
some  extent,  tunnel  across  the  Santa  Monica  mountains  to  enter  the  Los  Angeles 
basin.  Significant  energy  also  must  be  reradiated  as  body  wave  energy  when  the 
surface  waves  reach  the  far  side  of  the  basins  since  the  surface  waves  do  not 
either  continue  into  the  next  basin  or  reflect  back  from  the  edge.  Little  energy 
reverses  direction  and  travels  back  to  the  north  in  the  basins.  The  peak  ampli¬ 
tudes  of  the  velocity  traces  are  greater  for  the  realistic  structure  than  for  the 
flat-layered  structure  in  the  San  Fernando  basin,  but  the  reverse  is  true  in  the 
Santa  Monica  mountains.  In  the  Los  Angeles  basin,  the  peak  amplitudes  are  less 
than  in  the  flat-layered  case  because  the  receivers  are  in  the  shadow  of  the  moun¬ 
tains,  for  the  case  of  sources  S  and  I.  The  source  D,  however,  excites  the  Los 
Angeles  basin  with  a  direct  S  wave,  which  is  more  efficiently  converted  to  surface 
waves  in  the  realistic  than  the  flat-layered  geometry. 

Note  that  the  seismograms  for  the  source  D  in  Figure  6  show  many  of  the 
characteristics  of  the  data  displayed  in  Figure  2,  that  is,  reduced  amplitude  with 
simple  direct  pulses  at  the  Santa  Monica  mountains  and  significant  surface  wave 
arrivals  in  each  basin.  Source  I,  at  6  km  depth,  matches  the  data  nearly  as  well, 
but  has  too  much  energy  in  surface  waves  compared  to  direct  waves.  For  the 
rest  of  the  paper  we  will  consider  only  a  point  dislocation  source  at  10  km  depth. 
Most  probably,  the  source  was  actually  finite,  but  most  of  the  effect  of  source 
finiteness  may  be  included  in  the  source  time  function  for  the  purpose  of  under¬ 
standing  the  motions  in  the  basins.  What  we  do  not  include  seems  less 
significant  than  the  effect  of  the  structure,  which  is  the  focus  of  this  paper.  This 
conclusion  is  justified  by  noting  that  a  point  source  with  an  empirical  source  time 
function,  an  independent  velocity  structure,  and  a  teleseismically  determined 


mechanism  expiai.to  most  of  the  velocity  records  observed  in  the  basin.  This 
approach  introduces  a  slight  contradiction;  we  put  both  pulses  at  10  km  depth 
with  the  same  mechanism,  while  Langston  (1078)  and  Heaton  (1982)  separate  the 
two  pulses  and  give  them  distinct  mechanisms.  Either  the  separation  and 
dififering  mechanisms  are  less  than  were  estimated  or  our  simpler  approach  pro¬ 
duces  a  similar  result,  for  in  the  end  we  can  match  the  data  within  the  limita¬ 
tions  of  our  knowledge  of  the  structure  and  our  2-D  approach. 

The  FD  seismograms  in  Figure  7b  suggest  that  the  body  wave  appears  at  the 
ridge  with  a  time  function  very  similar  to  the  source  time  function.  The 
waveform  of  the  transverse  component  at  the  ridge  station  D068  is  shown  in  velo¬ 
city  and  displacement  in  Figure  8.  The  displacement  trace  shows  two  strong 
pulses  of  energy.  Heaton’s  (1982)  source  model  predicts  a  strong  second  arrival 
with  a  4  second  time  delay,  while  Langston  (1978)  finds  a  value  of  4.9  seconds. 
Both  models  appear  to  fit  the  observations  at  the  ridge  quite  well.  The  relative 
amplitude  ratio  of  the  two  pulses  appear  to  favor  the  Langston  source  model, 
which  is  given  in  Table  2. 

One  possible  strategy  at  this  point  would  be  to  add  the  secondary  source 
and  make  adjustments  in  the  source  time  description.  Another  approach  that 
might  be  called  an  "empirical  source  model"  is  to  assume  that  the  displacement 
record  D068  is  a  good  source  description  in  this  particular  direction  and  use  it  to 
predict  the  other  seismograms.  This  empirical  approach  is  adopted  and  the  FD 
seismograms  with  the  empirical  source  model  are  presented  in  Figure  9.  The 
moment,  which  is  difficult  to  control  because  the  long-period  information  is 
absent  from  the  data,  is  set  so  that  the  amplitude  of  the  observed  transverse 
record  at  station  D068  matches  the  synthetic  at  that  range.  When  a  cluster  of 
stations  are  at  nearly  the  same  range,  only  one  representative  trace  is  plotted. 
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In  Figure  9,  the  agreement  between  the  data  and  the  FD  seismograms  is 
excellent  at  all  but  the  closest  range,  C041  at  Pacoima  dam,  where  the  near-field 
assumptions  break  down  and  finiteness  and  directivity  may  become  important. 
The  frequency  content  of  the  records  in  the  San  Fernando  valley  is  higher  in  the 
data  than  in  the  synthetics,  which  is  probably  due  to  inaccuracies  in  the  velocity 
profile  which  produce  too  low  a  resonant  frequency  in  the  FD  modeling. 

The  P-SV  motions  can  be  simulated  by  constructing  Green’s  functions  as  for 
the  SH  case  except  that  two  rather  than  three  fundamental  fault  orientations  are 
required.  These  responses  for  the  assumed  two-dimensional  structure  and  10  km 
deep  source  are  given  in  Figure  10.  As  discussed  earlier,  the  vertical  radiation 
patterns  strongly  influence  the  relative  body-wave  to  surface-wave  ratios.  The 
A  i,  A  2,  and  A3  factors  are  all  significant  for  the  deep  source  mechanism.  The 
A  3  or  45°  dip-slip  component  dominates,  as  might  be  expected  for  a  mostly 
thrust  event,  because  A  3  is  slightly  larger  than  A  j  and  A  2  and  the  amplitudes  of 
the  Green’s  functions  for  A  3  are  larger  than  those  of  A  j  and  A  2. 

The  FD  seismograms  are  compared  with  the  data  for  the  vertical  and  radial 
components  in  Figures  11  and  12.  Again,  when  a  cluster  of  stations  are  at  nearly 
the  same  range,  only  one  representative  trace  is  plotted.  The  amplitude  of  the 
synthetic  seismograms  is  determined  from  the  SH  scaling  described  above.  The 
match  in  timing  is  somewhat  arbitrary  since  there  is  no  absolute  timing  for  the 
data. 

The  match  in  amplitude  between  the  P-SV  synthetic  velocities  and  the  data 
is  good.  In  addition,  the  duration  of  shaking  in  the  San  Fernando  basin  and  the 
waveforms  at  the  start  of  the  record  in  the  Los  Angeles  basin  are  also  well- 
modeled  using  the  empirical  time  function.  As  in  the  SH  case,  the  Santa  Monica 
mountains  produce  a  strong  decrease  in  the  amplitude  of  the  velocity  traces,  and 
both  basins  apparently  convert  the  direct  waves  into  surface  waves,  Rayleigh 


waves  in  this  ca™.  at  the  edge  nearest  the  source. 

Two  troublesome  problems  appear  in  the  forward  modeling  of  the  San  Fer¬ 
nando  records.  First,  the  small  amount  of  shallow  moment  release  allowed  by 
our  model  is  in  conflict  with  observations  of  1  meter  of  thrust  motion  measured 
at  the  surface  by  Alewine  (1974),  as  well  as  the  suggestions  of  Heaton  (1982), 
Langston  (1978)  and  Alewine  (1974)  of  significant  shallow  moment  release.  In  our 
model,  sources  at  between  6  and  14  km  depths  excite  basins  in  a  similar  way,  so 
it  is  only  in  the  top  5  km  that  we  would  exclude  significant  moment  release.  The 
half-space  model  of  Heaton  (1982)  and  the  layered  model  of  Langston  (1978) 
overestimate  the  velocity  near  the  surface,  so  a  smaller  moment  would  serve  to 
excite  the  same  energy  in  the  seismic  waves  given  the  correct  lower  velocities. 
Also,  it  is  possible  that  the  energy  radiated  from  the  shallow  portion  of  the  fault 
is  absorbed  by  some  mechanism  such  as  decoupling  of  the  two  sides  of  the  fault 
plane  or  very  low  Q  in  the  region,  but  these  explanations  are  not  compelling. 
The  main  point  of  this  paper  is  that  the  relative  amplitudes  and  duration  of 
shaking  across  basins  can  be  explained  with  FD  modeling,  and  the  details  of 
faulting  in  our  model  are  secondary.  To  the  extent  that  shallow  sources  tend  to 
excite  larger  motions  in  the  San  Fernando  basin  relative  to  those  in  the  Los 
Angeles  basin  than  were  observed,  this  study  suggests  less  moment  release  shal¬ 
lower  than  about  6  km  depth  than  the  studies  of  Heaton  (1982)  and  Langston 
(1978). 

The  second  problem  may  be  seen  in  Figures  11  and  12,  where  the  Rayleigh 
wave  created  in  the  San  Fernando  basin  tunnels  across  the  Santa  Monica  moun¬ 
tains  and  has  a  much  larger  amplitude  in  the  Los  Angeles  basin  than  is  seen  in 
the  data.  We  have  tried  rather  extreme  structures  such  as  more  separation 
between  the  basins  and  different  kinds  of  edges  on  the  basins,  but  the  Rayleigh 
waves  simply  are  better  at  tunneling  across  the  mountains  than  the  Love  waves, 


and  it  is  diffi  c  !■„  iiiatcli  the  data.  Three-dimensional  effects  may  be  responsi¬ 
ble  for  this  problem.  Since  this  profile  skirts  the  west  edge  of  the  Los  Angeles 
basin,  the  Rayleigh  waves  jumping  the  mountains  will  be  refracted  into  the 
slower  material  in  the  center  of  the  basin  and  the  amplitude  of  these  Rayleigh 
waves  at  stations  Hi  18  and  S267  may  be  small. 

An  alternative  explanation  may  be  that  the  velocity  model  for  the  San  Fer¬ 
nando  basin  has  too  thick  a  column  of  slow  sediments.  In  the  data,  it  appears 
that  the  surface  waves  are  shorter  period  in  the  San  Fernando  basin  than  the  Los 
Angeles  basin.  In  the  model,  the  same  3-5  second  period  surface  waves  are 
excited  in  both  the  San  Fernando  and  Los  Angeles  basins,  suggesting  that  the 
San  Fernando  basin  model  is  inaccurate  in  that  it  should  have  a  thinner  layer  of 
slow  sediments  than  the  Los  Angeles  basin.  If  the  two  basins  had  distinctly 
different  resonant  periods,  the  surface  waves  from  the  San  Fernando  basin  would 
not  excite  the  Los  Angeles  basin  as  much. 

The  peak  amplitude  comparison  between  the  data  and  the  FD  synthetics  is 
summarized  in  Figure  13.  The  match  is  very  good.  We  did  not  deviate  from  the 
velocity  model  derived  from  Duke  et  a!.  (1971)  to  keep  this  as  much  of  a 
forward-modeling  exercise  as  possible.  We  feel  as  a  forward-modeling  exercise, 
this  simulation  has  been  successful. 

For  insight  into  the  creation  and  destruction  of  the  surface  waves  at  the 
edges  of  the  basins,  the  envelope  of  the  transverse  velocity  records  is  shown  in 
Figure  14.  For  this  plot,  250  rather  than  eight  receivers  are  used.  In  this  figure, 
a  sketch  of  the  basins  and  mountains  is  included  at  the  bottom  for  location.  We 
use  the  same  source  location  and  empirical  time  function  as  we  did  for  the  previ¬ 
ous  synthetics  seismograms.  The  first  arrival  in  the  seismograms  displayed  in 
Figures  9  and  14  is  the  direct  SH  wave.  The  direct  wave  advances  slightly  at  the 
Santa  Monica  mountains  because  of  the  faster  material  at  the  surface.  In  both 


the  San  Fernando  basin  and  the  Los  Angeles  basin,  the  direct  SH  wave  incident 
upon  the  edge  of  the  basin  nearest  the  epicenter  produces  a  surface  wave  train 
(Love  waves)  that  crosses  the  basin  then  converts  back  to  body  waves. 

Some  energy  from  the  Love  wave  in  the  San  Fernando  basin  converts  to  a 
body  wave  with  a  rapid  apparent  velocity  across  the  Santa  Monica  mountains, 
then  partially  converts  back  to  a  surface  wave  in  the  Los  Angeles  basin.  This 
arrival  may  be  thought  of  as  surface  wave  energy  that  has  tunneled  across  the 
mountains.  The  surface  wave  generated  by  the  direct  SH  pulse,  however,  is 
stronger  than  the  tunneling  surface  wave  for  this  source  location  and  mechanism. 

CONCLUSIONS 

The  strong  motions  recorded  within  the  San  Fernando  and  Los  Angeles 
basins  during  the  San  Fernando  earthquake  of  1971  can  be  largely  understood  by 
forward  modeling  with  teleseismically  determined  source  parameters  through 
known  structure  with  only  two  dimensional  variations  in  velocity  and  density. 
This  modeling  will  accept  less  moment  release  shallower  than  5  km  depth  than 
has  been  suggested  in  some  previous  studies,  but  this  inconsistancy  might  be 
reconciled  if  the  fully  three  dimensional  structure  appropriate  for  closed  basins  is 
used  in  the  FD  calculation.  The  dramatic  differences  between  using  the  laterally 
varying  structure  of  Duke  et  al.  (1971)  and  a  flat-layered  structure  render  it 
impossible  to  fit  this  data  from  the  San  Fernando  earthquake  with  a  flat-layered 
model. 

The  geometry  of  the  basin  structure  is  important  in  determining  the 
attenuation  of  peak  velocity  along  a  profile,  as  has  been  previously  suggested  by 
theoretical  studies  (Boore,  1970,  Lysmer  and  Drake,  1972,  Bard  and  Gariel,  1986). 
The  mountains  shadow  the  basins  behind  them  against  shallow  surface  waves, 
but  the  edges  of  the  basin  nearest  the  earthquake  tends  to  convert  body  waves 
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again  to  shallc"-  surface  waves.  The  source  depth  of  the  earthquake  is  less 
important  than  for  the  flat-layered  case  in  determining  the  duration  of  shaking 
and  attenuation  of  peak  amplitude  with  distance.  Basins  with  sharp  edges  tend 
to  generate  some  back-scattered  surface  waves;  basins  with  gradual  edges  have 
negligible  back-scattered  energy. 


It  is  encouraging  that  the  2-10  second  period  seismic  waves  observed  in  the 
San  Fernando  and  Los  Angeles  basins  can  be  largely  explained  by  known  struc¬ 
ture.  We  feel  that  this  line  of  deterministic  modeling  can  help  us  understand  and 
reduce  the  uncertainties  in  amplitude  attenuation  with  distance  of  strong 
motions. 
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Table  1  :  Velocity  and  density  structure  for  the  cross-section  in  the  middle  of  the  San 
Fernando  basin 


P-wave 

S-wave 

Density 

velocity 

velocity 

1.2 

0.6 

1.7 

2.0 

1.1 

1.8 

2.5 

1.4 

1.9 

3.1 


1.8 


2.1 


Table  2  :  Source  parameters  of  the  San  Fernando  earthquake 


Study 

Method 

Strike 

Dip 

Rake 

Depth 

km 

Whitcombf  19731 

-67° 

52° 

72° 

8 

First  motions 

Alewinef 19741 

-67° 

53° 

72° 

Static  displacements 

Alewinef  19741 

-67° 

53° 

66-82° 

Surface  waves 

Lanyston(1978) 

first  source 

-79° 

44° 

15 

Teleseismic  records 

second  source 

-80° 

18° 

96° 

Lanestonf  19781 

first  source 

-70° 

53° 

76° 

8-15 

Teleseismic  records 

second  source 

-80° 

29° 

90° 

Heatonfl982l  first  source 


combined  study  second  source 
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FIGURE  CAPTIONS 

Figure  1.  Map  and  cross-section  of  the  San  Fernando  region  from  Duke  et 
al.  (1971).  The  epicenter  is  marked  by  a  cross.  The  surface  breakage  is 
indicated  by  the  cross-hatched  line.  The  filled  triangles  are  the  locations 
of  the  strong-motion  instruments  used  in  this  section.  Cross-hatched  areas 
show  surface  exposure  of  bedrock.  The  bottom  of  the  basin  for  the  profile 
A-A’  is  shown  below,  where  dashed  portions  show  where  the  boundary  is 
not  known.  The  cross-section  has  vertical  exageration  of  2:1.  SFB  indi¬ 
cates  the  San  Fernando  basin,  SMM  the  Santa  Monica  mountains,  LAB 
the  Los  Angeles  basin,  PVH  the  Palos  Verdes  hills,  and  PO  the  Pacific 
ocean. 

Figure  2a.  Velocity  records  of  the  1971  San  Fernando  earthquake,  taken 
from  EERL  (1974).  The  traces  are  aligned  relative  to  a  high-frequency, 
early  arrival  on  the  vertical  component  that  is  interpreted  to  be  a  direct 
compressional  wave.  Amplitude  is  given  in  cm/sec.  The  station  names 
are  listed  at  the  far  left.  The  stations  are  shown  in  order  of  increasing  epi- 
central  distance,  but  the  actual  station  spacing  is  irregular.  Figure  2b. 
Smoothed  velocity  records.  The  records  shown  in  Figure  2a  are  convolved 
by  a  gaussian  pulse  about  one  second  wide  to  filter  out  frequencies  that 
cannot  be  properly  handled  by  the  finite-difference  grid. 

Figure  3.  Stratigraphic  and  structural  cross-section  of  the  profile  from  the 
San  Gabriel  mountains  on  the  left,  across  the  San  Fernando  basin,  the 
Santa  Monica  mountains,  and  the  Los  Angeles  basin,  to  the  Palos  Verdes 
hills  on  the  right.  The  stratigraphy,  the  velocities,  and  the  densities  are 
taken  from  Duke  et  al.  (1971),  except  that  the  velocities  below  the  basins 
are  taken  from  Kanamori  and  Hadley  (1975).  The  heavy  vertical  lines  in 
the  stratigraphy  diagram  indicate  the  location  of  well  logs  used  by  Duke  et 
al.  (1971)  to  construct  the  model.  Only  the  shear  wave  velocity  is  shown, 
the  compressional  wave  velocity  and  the  density  corresponding  to  each 
shear  wave  velocity  are  given  in  Table  1.  The  letters  B,  D,  I,  and  S 
correspond  to  the  locations  of  the  bottommost,  deep,  intermediate,  and 
shallow  sources  discussed  in  the  text.  The  triangles  along  the  surface 
mark  the  location  of  the  receivers  for  the  finite  difference  seismograms. 
The  vertical  exaggeration  is  2:1. 

Figure  4.  Schematic  north-south  cross-section  of  results  from  teleseismic 
body  wave  studies.  Langston  (1978)  interprets  the  faulting  of  the  San  Fer¬ 
nando  earthquake  as  two  point  sources  in  Figure  4a  and  as  two  finite  seg¬ 
ments  in  Figure  4b.  The  moment  is  assumed  to  be  uniformly  distributed 
with  depth  for  each  source  in  Figure  4b.  Heaton  (1982)  combines  the 
teleseismic  data,  near-field  data  on  hard  rock  sites,  and  static  data  to  find 
the  two  non-intersecting  finite  fault  segments  shown  in  Figure  4c.  The 
moment  is  assumed  to  be  non-uniformly  distributed  with  depth  for  each 


source  in  Figure  4c.  Details  of  these  models  may  be  found  in  Table  2. 

Figure  5.  Green’s  functions  for  the  transverse  component  of  velocity  with 
S  and  I  source  locations.  The  strike-slip  case  corresponds  to  _A4=  1  and 
>t6  =  o,  and  the  dip-slip  case  corresponds  to  Ai  =  0  and  At  —  l.  Moment 
scaling  is  discussed  in  the  text.  The  source  time  function  is  impulsive,  and 
the  frequency  content  is  limited  by  the  seismic  Q  of  25. 

Figure  6.  Green’s  functions  for  the  transverse  component  of  velocity  with 
the  D  and  B  source  locations.  The  strike-slip  case  corresponds  to  A  4  -=  1 
and  j46  =  0,  and  the  dip-slip  case  corresponds  to  A4  =  0  and  ^46=1. 
Moment  scaling  is  discussed  in  the  text.  The  source  time  function  is 
impulsive,  and  the  frequency  content  is  limited  by  the  seismic  Q  of  25. 

Figure  7a.  Finite-difference  seismograms  for  the  transverse  component  of 
velocity  for  the  model  with  flat  layers  of  the  thickness  and  velocities 
appropriate  for  the  middle  of  the  San  Fernando  basin  in  Figure  3,  listed  in 
Table  1.  The  source  has  a  gaussian  time  function  3.0  seconds  in  width 
that  is  shown  in  Figure  8c.  The  mechanism  is  strike-slip  and  the  moment 
scaling  is  discussed  in  the  text.  Velocities  are  in  cm/sec.  Figure  7b.  Fin¬ 
ite  difference  velocity  seismograms  for  the  model  shown  in  Figure  3.  The 
source  has  the  same  time  function  3.0  seconds  in  width  as  Figure  7a.  The 
mechanism  is  strike-slip.  The  location  of  the  shallow,  intermediate,  and 
deep  depth  sources  are  indicated  in  Figure  3. 

Figure  8.  The  calculation  of  an  empirical  source  time  function  for  the 
deep  source.  Figure  8a  shows  the  observed  transverse  velocity  waveform 
at  the  Santa  Monica  mountains  station  D068.  Figure  8b  shows  the  dis¬ 
placement  waveform  at  station  D068.  Predictions  of  secondary  arrivals  by 
Heaton  (1982)  and  Langston  (1978)  are  indicated  by  the  H  and  L  above 
the  trace.  Figure  8c  shows  the  gaussian  time  function  assumed  in  the  con¬ 
struction  of  the  synthetics  shown  in  Figure  7. 

Figure  9.  Comparison  of  filtered  transverse  component  of  data  with  FD 
seismograms  computed  with  the  empirical  time  function  shown  in  Figure 
8.  The  heavy  traces  show  the  data,  with  the  station  name  to  the  left  and 
the  amplitude  in  cm/sec  to  the  right.  The  light  traces  show  the  FD 
seismograms,  with  the  range  in  km  to  the  left  and  the  amplitude  to  the 
right.  The  FD  seismograms  assume  point  source  D  in  the  structure  shown 
in  Figure  3,  and  a  Q  of  25.  The  source  time  function  is  chosen  and  the 
source  amplitude  is  scaled  so  that  the  waveform  and  amplitude  of  the  syn¬ 
thetic  seismogram  at  the  ridge  station  (30  km)  matches  that  of  station 
D068. 

Figure  10.  Vertical  and  radial  Green’s  functions  for  source  D.  The 
strike-slip  case  corresponds  to  ^  i  «•  l,  /42«=0,  and  — 0.  The  dip-slip 
case  corresponds  to  A,*0,  A3=i,  and  y4,  =  0.  The  45*  dip-slip  case 
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corresponds  to  A ,  =  0,  At  =  0,  and  As  =  1.  The  moment  for  each  Green’s 
function  is  lO2®  dyne-cm.  The  frequency  content  is  limited  mostly  by  the 
seismic  Q  of  25. 

Figure  11.  Comparison  of  filtered  vertical  component  of  data  with  FD 
seismograms  computed  with  the  empirical  time  function  shown  in  Figure 
8.  The  heavy  traces  show  the  data,  with  the  station  name  to  the  left  and 
the  amplitude  in  cm/sec  to  the  right.  The  light  traces  show  the  FD 
seismograms,  with  the  range  in  km  to  the  left  and  the  amplitude  to  the 
right.  The  data  is  lined  up  on  a  high-frequency  arrival  on  the  vertical 
component  that  is  presumably  a  direct  compressional  wave.  The  FD 
seismograms  assume  a  point  source  at  the  location  D  (10  km  depth)  in  Fig¬ 
ure  3;  they  are  computed  with  the  structure  shown  in  Figure  3,  and  a  Q  of 
25  is  assumed.  The  amplitudes  of  the  finite  difference  seismograms  are 
consistent  with  those  for  the  transverse  component  in  Figure  9. 

Figure  12.  Comparison  of  filtered  radial  component  of  the  data  with  FD 
seismograms  computed  with  the  empirical  time  function  shown  in  Figure 
8.  The  heavy  traces  show  the  data,  with  the  station  name  to  the  left  and 
the  amplitude  in  cm/sec  to  the  right.  The  light  traces  show  the  FD 
seismograms,  with  the  range  in  km  to  the  left  and  the  amplitude  to  the 
right.  The  data  is  lined  up  on  a  high-frequency  arrival  on  the  vertical 
component  that  is  presumably  a  direct  compressional  wave.  The  FD 
seismograms  assume  a  point  source  at  the  location  D  (10  km  depth)  in  Fig¬ 
ure  3;  they  are  computed  with  the  structure  shown  in  Figure  3,  and  a  Q  of 
25  is  assumed.  The  amplitudes  of  the  finite  difference  seismograms  are 
consistent  with  those  for  the  transverse  component  in  Figure  9. 

Figure  13.  Peak  velocity  attenuation  with  distance.  The  crosses  show  the 
peak  velocity  of  the  smoothed  data  versus  offset.  The  line  shows  the 
attenuation  for  the  FD  simulation  of  the  motions. 


Figure  14.  A  seismic  section  of  the  envelope  of  the  transverse  component 
of  velocity  across  the  San  Fernando  and  Los  Angeles  basins.  Dark  por¬ 
tions  of  the  images  have  energy  in  the  velocity  traces;  light  portions  have 
little  or  no  energy  at  that  time.  250  receivers  are  used,  rather  than  the 
eight  shown  up  to  this  point.  The  same  source  is  used  as  in  Figure  8.  The 
major  geologic  structures  are  sketched  below  the  section,  where  SFB  indi¬ 
cates  the  San  Fernando  basin,  SMM  the  Santa  Monica  mountains,  and 
LAB  the  Los  Angeles  basin. 
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Abstract 


The  methods  for  Representation  Theorem  (RT)  coupling  of  finite  element  (FE) 
or  finite  difference  calculations  and  Harkrider’s  (Harkrider  1964,  1970)  propagator 
matrix  method  calculations  to  produce  a  hybrid  method  for  propagation  of  SH  mode 
sum  seismograms  across  paths  that  contain  regions  of  non  plane-layered  structure  are 
explained  and  developed.  The  coupling  methods  explained  in  detail  use  a  2-D  Carte¬ 
sian  FE  formulation.  Analogous  methods  for  the  3-D  method  follow  directly.  Exten¬ 
sive  tests  illustrating  the  validity  and  accuracy  of  the  implementation  of  these  cou¬ 
pling  methods  are  discussed.  These  hybrid  techniques  are  developed  to  study  the 
propagation  of  surface  waves  across  regional  transition  zones  or  other  heterogeneities 
that  exist  in  part  of  a  longer,  mostly  plane-layered,  path.  The  effects  of  a  thinning  or 
thickening  of  the  crustal  layer  on  the  propagation  of  Lg  mode  sum  seismograms  have 
been  examined  in  this  study.  The  thinning  or  thickening  of  the  crustal  layer  is  used 
as  a  simple  model  of  ocean  continent  transitions.  The  Lg  phase  is  of  particular 
interest  since  it  is  used  in  several  important  applications  such  as  mapping  the  extent 
of  continental  crust,  magnitude  determination,  and  discrimination  between  explosive 
and  earthquake  sources.  The  understanding  of  the  observations  that  Lg  wave  is 
attenuated  completely  when  the  propagation  path  includes  an  oceanic  portion  of 
length  greater  than  one  hundred  to  two  hundred  kilometers  or  a  region  of  complex 
crustal  structure  is  not  complete,  and  a  clear  explanation  of  these  phenomena  could 
have  important  consequences  for  all  these  types  of  studies.  The  transition  model  cal¬ 
culations  done  in  this  study  show  that  passage  through  a  region  of  thinning  crustal 
thickness,  the  model  for  a  continent  to  ocean  transition,  increases  the  amplitude  and 
coda  length  of  the  Lg  wave  at  the  surface,  and  allows  much  of  the  modal  energy 


r.^rj  WVFTO  KF H"  JL*  Wit"  ^  *" 


-46- 


trapped  is  the  crust,  which  forms  the  Lg  phase,  to  escape  into  the  subcrustal  layers 
as  body  waves  or  other  downgoing  phases.  The  magnitude  of  both  these  effects 
increases  as  the  length  of  the  transition  increases  or  the  slope  of  the  layer  boundaries 
decrease.  The  passage  of  the  wavefront  exiting  the  continent  to  ocean  transition 
region  through  the  oceanic  structure  allows  further  energy  to  escape  from  the  crustal 
layer,  and  produces  a  decrease  in  Lg  amplitude  at  the  surface  as  the  length  of  the  oce¬ 
anic  path  increases.  The  amplitude  decrease  is  maximum  near  the  transition  region 
and  decreases  with  distance  from  it.  Passage  through  a  region  of  thickening  crust, 
the  model  of  a  ocean  to  continent  transition,  causes  a  rapid  decrease  in  the  Lg  ampli¬ 
tude  at  the  surface  of  the  crust.  The  energy  previously  trapped  in  the  oceanic  crustal 
layer  spreads  throughout  the  thickening  crustal  layer,  and  any  amplitude  which  has 
been  traveling  through  the  subcrustal  layer  but  has  not  reached  depths  below  the 
base  of  the  continental  crust  is  transmitted  back  into  the  continental  crust.  The 
attenuation  of  Lg  at  the  crustal  surface  along  a  partially  oceanic  path  occurs  in  the 
oceanic  structure  and  in  the  ocean  to  continent  transition  region.  The  attenuation 
at  the  surface  depends  in  part  on  the  escape  of  energy  at  depth  through  the  continent 
to  ocean  transition  region  into  the  underlaying  half-space.  The  total  attenuation  of 
Lg  due  to  propagation  through  a  forward  transition  followed  by  a  reverse  transition 
is  at  most  a  factor  of  four  to  six.  This  is  inadequate  to  explain  the  observed  attenua¬ 
tion  of  Lg.  Thus,  additional  effects,  other  than  geometry  must  be  considered  to  pro¬ 
vide  a  complete  explanation  of  the  attenuation  of  Lg. 
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Introduction 

This  thesis  presents  a  study  of  the  propagation  of  L(  waves  across  ocean  con¬ 
tinent  transition  regions.  The  transition  regions  are  represented  by  a  simplified  model 
consisting  of  a  crustal  layer  with  thickness  that  increases  or  decreases  smoothly 
between  oceanic  and  continental  values.  The  distance  over  which  this  smooth  varia¬ 
tion  occurs  is  referred  to  as  the  length  of  the  transition  region.  The  changes  to  a  L, 
wavefield,  when  it  travels  through  one  of  these  simple  transition  regions,  are  modeled 
using  a  hybrid  method  which  is  explained  and  developed  in  sections  four,  five,  and  six 
of  chapter  one.  The  hybrid  method  combines  the  Finite  Element  method  and  the  Pro¬ 
pagator  Matrix  technique  through  application  of  the  Seismic  Representation  Theorem. 
These  two  methods  are  explained  and  the  notation  associated  with  them  is  defined  in 
section  two  and  three  of  chapter  one.  In  the  regions  of  complex  structure,  those 
regions  containing  the  transition  regions,  the  Finite  Element  method  is  used  to 
transmit  the  Lg  waves.  In  the  regions  of  plane  layered  structure  either  the  Propagator 
Matrix  technique  is  used  to  directly  determine  the  seismograms  for  a  source  in  the 
layers  being  considered,  or  it  is  used  to  determine  the  Green’s  functions  for  propaga¬ 
tion  from  the  interface  between  the  complex  region  and  the  layered  region  to  the 
receiver.  In  the  first  case,  displacement  time  histories  are  determined  at  a  set  of  points 
equally  spaced  in  depth  along  a  vertical  plane.  These  time  histories  are  used  as  boun¬ 
dary  conditions  for  the  Finite  Element  method.  Extensive  tests  of  the  accuracy  of 
producing  a  Finite  Element  wavefield  by  specifying  the  displacements  as  a  function  of 
time  along  the  vertical  plane  defining  the  Finite  Element  grid  edge  are  presented  in 
sections  seven  and  eight  of  chapter  one.  The  results  of  these  tests  are  useful  in  under¬ 
standing  the  sources  of  uncertainty  due  to  truncation  of  the  vertical  plane  at  some 
depth  less  than  infinity.  In  the  case  where  the  propagator  method  is  used  to  deter¬ 
mine  Green’s  functions  the  Green’s  functions  are  combined  with  the  stress  and  dis¬ 
placement  time  histories  recorded  at  the  boundary  between  the  plane  layered  region 


and  the  complex  region  according  to  the  Representation  Theorem.  Tests  of  the  vali¬ 
dity  and  accuracy  of  this  method  of  coupling  are  discussed  in  section  two  of  chapter 
three.  The  results  of  the  study  on  the  effects  of  transition  length  on  the  propagation 
of  Lg  waves  are  reported,  discussed  and  interpreted  in  chapter  two.  The  results  of  the 
study  of  the  effects  of  the  length  of  the  intermediate  oceanic  path  on  the  propagation 
of  L(  are  mentioned  in  chapter  two  and  again  in  the  first  and  last  sections  of  chapter 
three.  In  order  to  understand  why  the  models  used  were  chosen  for  this  study,  and 
why  this  study  is  important,  it  is  useful  to  discuss  and  review  the  the  different  types 
of  studies  that  use  Lg  data.  In  the  remainder  of  this  section  the  background  of  the 
problem  will  be  discussed. 

In  order  to  choose  a  reasonable  form  for  a  simple  model  of  a  ocean  continent  or 
continent  ocean  transition  region  it  is  useful  to  examine  existing  observations  of  such 
regions.  Structures  of  transition  regions  have  been  studied  using  many  techniques, 
both  geological  and  geophysical.  These  techniques  include  reflection  and  refraction 
seismic  surveys,  gravity  studies,  and  geologic  mapping  and  interpretation.  Many  such 
studies  have  been  completed  so  only  a  few  examples  will  be  cited  here  to  support  the 
choice  of  the  model  structure  used.  A  review  of  studies  of  the  continental  margins  of 
the  eastern  and  western  coasts  of  Canada  by  Keen  and  Hyndman  (1979)  shows  sample 
sections  of  the  transition  regions  along  the  passive  margins  of  the  eastern  coast.  The 
transition  lengths  observed  vary  between  about  fifty  and  two  hundred  kilometers  with 
the  average  being  about  one  hundred  kilometers.  Sheridan  et  al.  (1979)  and  Grow  et 
al.  (1979)  discuss  refaction  studies  off  the  northeastern  United  States  coast  which 
show  similar  transition  geometries  and  transition  lengths  of  one  to  three  hundred 
kilometers.  Uchupi  and  Austin  (1979)  examine  reflection  profiles  off  the  eastern  coast 
of  North  America  which  indicate  that  the  length  of  the  upper  surface  of  the  transition 
region  is  about  one  hundred  kilometers.  Hinz  et  al.  (1979)*present  the  same  types  of 
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results  off  Labrador  and  Greenland  showing  transition  lengths  of  about  one  hundred 
kilometers.  Studies  off  northern  Alaska  by  Eittreim  and  Grantz  (1979)  and  Grantz  et 
al.  (1979)  show  transition  lengths  of  the  upper  transition  surface  of  between  fifty  and 
one  hundred  kilometers.  Barton  et  al.  (1984)  derive  a  structure  beneath  the  North  Sea 
which  shows  similar  transition  regions.  Le  Douaran  et  al.  (1984)  examine  the  crustal 
structure  in  the  northwestern  Mediterranean  Basin,  and  Makris  et  al.  (1983)  and 
Makris  and  Stobbe  (1984)  derive  structures  in  the  Eastern  Mediterranean.  In  the 
Mediterranean  region  the  transition  regions  are  fifty  to  one  hundred  fifty  kilometers  in 
length.  Most  of  the  transition  region  structures  mentioned  above  show  similarities.  In 
all  cases  the  crustal  layer  thins  from  the  continent  toward  the  ocean.  In  most  cases 
this  thinning  is  fairly  constant.  The  thicknesses  of  the  continental  and  oceanic  crustal 
layers  vary  somewhat.  The  resolution  and  detail  in  the  models  vary  considerably.  The 
details  would  need  to  be  accounted  for  if  cases  of  detailed  waveform  modeling  were  to 
be  undertaken,  but,  they  are  not  critical  for  this  study.  The  lengths  of  transition 
regions  suggested  by  the  observations  are  between  fifty  and  three  hundred  kilometers. 
These  studies  all  indicate  that  the  simplified  models  used  in  this  investigation  are  rea¬ 
sonable.  These  models  use  a  thirty  two  kilometer  crustal  layer  over  a  half-space  as  a 
model  of  the  continental  structure  and  a  five  kilometer  ocean  over  a  five  kilometer 
crust  over  a  half-space  as  a  model  of  the  oceanic  structure.  The  thickness  of  the  crus¬ 
tal  layer  varies  smoothly  between  the  oceanic  and  continental  values  as  one  passes 
through  a  transition  region.  Model  transition  regions  are  constrained  by  computa¬ 
tional  limits  to  lengths  of  one  hundred  kilometers  or  less. 

The  wavefields  transmitted  through  these  models  are  synthesized  to  model  tLe  L, 
arrivals  in  a  seismogram.  The  Lg  arrival  was  first  defined  by  Press  and  Ewing  (1952) 
Based  on  studies  of  seismograms  recorded  in  North  America,  they  defined  L#  to  have 
periods  of  0.5  to  6  seconds,  phase  velocities  between  2.0  and  3.5  km/s,  reverse 


dispersion  at  distances  greater  than  20°,  and  a  sharp  high  amplitude  onset.  They 
noted  that  Lg  arrivals  are  visible  for  purely  continental  paths  and  are  gradually  elim¬ 
inated  as  the  oceanic  portion  of  the  path  increases  in  length  beyond  one  hundred 
kilometers.  The  arrivals  are  visible  on  all  three  components  but  are  larger  on  the  hor¬ 
izontals.  They  suggested  two  possible  interpretations  of  Lg,  then  abandoned  the  one 
which  explained  Lt  in  terms  of  higher  mode  Love  waves  in  favor  of  one  which 
explained  them  as  channel  waves  trapped  in  the  crust  above  a  low  velocity  layer. 
The  idea  that  Lg  was  a  crustal  channel  wave  was  expanded  and  discussed  by  Press 
and  Gutenburg  (1956)  and  Gutenburg  (1955)  but  even  they  stated  that  such  an  explsr 
nation  was  at  best  tentative.  The  interpretation  of  Lg  in  terms  of  higher  mode  Love 
waves  was  abandoned  by  Press  and  Ewing  (1952)  because  of  the  motions  on  the  verti¬ 
cal  and  longitudinal  components  and  the  long  coda.  Oliver  and  Ewing  (1957,  1958a, 
1958b)  and  Oliver  et  al.  (1959)  studied  higher  mode  Love  waves  and  Rayleigh  waves 
and  concluded  that  the  vertical  and  longitudinal  components  of  Lg  are  composed  of 
higher  mode  Rayleigh  waves  and  that  higher  mode  Love  waves  could  account  for  the 
transverse  component.  They  considered  only  the  first  two  higher  modes  and  thus  did 
not  explain  the  long  coda,  but  only  its  initial  portion.  Kovach  and  Anderson  (1964) 
showed  in  detail  that  Lg  could  be  interpreted  using  higher  mode  Rayleigh  and  Love 
waves.  Panza  et  al.  (1972)  related  the  higher  mode  Rayleigh  wave  interpretation  to 
the  channel  wave  explanation.  Panza  and  Calcagnile  (1974,  1975)  compare  models 
with  and  without  a  low  velocity  channel  below  the  crustal  waveguide  and  conclude 
that  such  a  channel  is  not  necessary  for  the  existence  of  Lg  waves.  Thus,  the  Lg 
arrival  is  clearly  explained  as  a  superposition  of  the  higher  mode  Love  waves  and  Ray¬ 
leigh  waves. 

All  the  studies  of  higher  mode  Rayleigh  waves  and  Love  waves  discussed  above 
are  based  on  analysis  of  group  velocities,  and  phase  velocities,  with  the  Lg  arrivals 


corresponding  to  the  group  velocities  of  the  extrema  of  the  group  velocity  curves. 
Knopoff  et  al.  (1975)  used  the  higher  mode  interpretation  of  Lg  to  calculate  synthetic 
seismograms  for  Lg.  These  seismograms  further  establish  the  validity  of  the  previous 
interpretations  based  on  group  velocities  of  higher  mode  Love  and  Rayleigh  waves  by 
matching  not  only  arrival  times  but  also  relative  amplitudes.  These  seismograms  also 
establish  that  the  observed  properties  of  L,  waves  can  be  produced  without'including 
a  low  velocity  channel.  Knopoff  et  al.  (1975)  continued  their  studies  using  synthetic 
L,  waves  to  define  properties  of  Lg  that  can  be  used  to  discriminate  differences  in  crus¬ 
tal  structure.  Bouchon  (1981,1982)  used  an  alternate  method  for  determining  Lg  syn¬ 
thetics  at  short  distances  (150  •  350  km)  to  model  data  from  an  earthquake  in  France. 
He  concluded  that  the  Lg  arrivals  were  composed  of  multiply  reflected  post-critical  SH 
and  SV  rays  and  that  synthetics  modeled  the  data  well  for  group  velocities  between 
3.5  and  2.8  km/s.  Several  other  workers  have  used  the  fact  that  Love  waves  can  be 
expressed  as  the  superposition  of  post-critical  multiply  reflected  SH  rays  in  the  crust 
to  investigate  the  Lg  arrivals.  Herrin  and  Richmond  (1960)  along  with  corrections  in 


Herrin  (1961)  used  ray  theory  to  show  Lg  could  be  modeled  using  multiply  reflected 
and  refracted  ray  in  a  thin  crustal  layer.  This  type  of  approach  was  also  used  by  Pec 
(1967)  and  Kennett  (1986)  to  address  some  properties  of  Lg.  Cara  and  Minster  (1981) 
developed  a  method  to  analyze  Rayleigh  type  L,  seismograms  recorded  along  a  linear 
array.  Cara  et  al.  (1981)  applied  this  method  to  data  recorded  in  the  northwestern 
Sierra  Nevada  and  in  southern  California.  They  conclude  that  representation  of  Lg  as 
a  single  multimode  wavetrain  is  an  oversimplified  approach  which  can  explain  observa¬ 
tions  in  the  initial  portions  of  the  wavetrain  (group  velocity  >  3.2  km/s),  while  the 
later  portions  of  the  waveform  are  strongly  affected  by  reflections  and  diffractions  from 
any  discontinuities  in  structure.  Other  phases  previously  defined  in  terms  of  the  chan¬ 
nel  model  have  also  been  modeled  using  the  higher  mode  surface  wave  model.  Schwab 
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et  &1.  (1974),  Mantovani  et  al.  (1977),  Nakanishi  et  al.  (1977),  and  Mantovani  (1978) 
considered  SB,  Stephens  and  Isacks  (1977)  considered  the  transverse  component  of  SB, 
and  Cansi  and  Bethoux  (1985)  investigate  the  T  phase  which  converted  to  Lg  at  the 
continental  margin.  In  light  of  all  the  evidence  the  higher  mode  interpretation  of  Lg 
best  explains  the  observed  properties  of  Lg,  even  though  variations  in  the  crustal 
structure  may  cause  a  longer  coda  than  that  calculated  using  a  simple  layered  struc¬ 
ture.  The  SH  Lg  seismograms  used  to  produce  the  wavefields  used  to  drive  the  FE  cal¬ 
culations  in  this  study  are  determined  as  mode  sums  over  the  fundamental  and  first 
five  higher  modes. 

One  of  the  important  uses  of  Lg  has  been  to  distinguish  regions  with  oceanic 
crustal  structures  from  those  with  continental  crustal  structures  and  to  map  regions 
where  sudden  changes  in  crustal  structure  occur.  While  the  theoretical  interpretation 
of  Lg  were  being  debated  many  observational  studies  of  Lg  waves  were  being  con¬ 
ducted.  These  studies  were  conducted  to  map  regions  of  oceanic  crustal  structure  and 
regions  of  complex  crustal  structure  on  continents  which  caused  the  Lg  arrival  to  be 
heavily  attenuated  or  vanish.  The  studies  considered  a  strong  Lg  to  indicate  con¬ 
tinental  structure  and  a  missing  or  highly  attenuated  Lg  to  indicate  passage  through 
oceanic  crust.  Intermediate  cases  were  included  in  one  group  or  the  other  depending 
on  the  workers,  but  none  used  a  quantitative  scale. 

Bath  (1954,1958)  collected  a  large  sample  of  Lg  data  along  Euroasiatic  paths  and 
fashioned  explanations  of  his  data  based  on  the  guided  wave  interpretation  of  Lg.  He 
also  defined  two  arrivals  be  called  Lg|  and  L^  with  distinctly  different  phase  velocities 
of  3.54  km/s  and  3.37  km/s.  His  observations  supported  the  conclusion  of  Press  and 
Ewing  (1952)  that  partially  oceanic  paths  resulted  in  attenuation  of  the  Lg  phase. 
However,  he  interpreted  some  portions  of  the  Arctic  Ocean  where  Lg  was  not 
attenuated  as  being  continental  in  structure.  This  approach,  of  defining  paths  which 


pass  under  oceans  but  do  not  attenuate  L(  as  continental  paths  has  also  been  used  in 
most  of  the  more  recent  studies  of  this  type  since  that  time.  These  studies  used  Lg 
recorded  in  many  regions  of  the  world.  Press  et  al.  (1956)  considered  paths  in  Africa, 
Press  (1956)  used  paths  in  California,  Savarensky  and  Valdner  (1960)  studied  the 
Black  Sea  region,  Bolt  (1957)  studied  data  from  Australia,  Utsu  (..)  considered  Japan 
and  the  region  surrounding  it,  and  Lehmann  (1952,  1957)  studied  significant  events  in 
North  America  and  Europe.  Oliver  et  al.  (1955)  considered  L,  propagation  in  the  Arc¬ 
tic  region,  Herrin  and  Minton  (1960)  in  the  Southwestern  United  States  and  Mexico, 
and  Wetmiller  (1974)  investigated  the  crustal  structure  in  the  Baffin  Bay  area.  Gre- 
gersen  (1984)  studied  the  crustal  structure  near  Denmark  and  in  the  North  Sea.  Addi¬ 
tional  data  for  this  area  were  presented  by  Kennett  and  Mykkeltveit  (1984)  who  used 
the  supplemented  data  set  and  the  method  of  Kennett  (1984a)  to  produce  synthetics 
and  an  interpretation  of  the  data  in  terms  of  a  multimode  wavetrain  traveling  through 
a  crustal  ’pinch’. 

As  originally  noted  by  Bath,  some  purely  continental  paths  show  attenuated  L( 
and  the  later  studies  mapped  these  regions.  Gumper  and  Pomeroy  (1970)  made  a 
more  extensive  study  of  L,  and  other  phases  in  Africa.  They  found  L(  velocities 
between  3.48  and  3.60  km/s  and  a  discontinuity  in  the  crust  near  the  rift  zone. 
Ruzaikin  et  al.  (1977)  mapped  regions  of  L#  propagation  in  Asia  and  saw  attenuation 
in  Tibet  and  regions  of  China.  Kadinsky-Cade  et  al.  (1981)  conducted  a  detailed 
study  of  L,,  S„  Pg,  and  PB  propagation  in  India,  the  Himalayan  Arc  and  the  Tibetan 
Plateau  producing  a  more  detailed  map  than  Ruzaikin  et  al.,  and  agreeing  with  their 
observations  and  the  sparse  observations  of  Bath  that  showed  the  absence  of  L,  on 
paths  crossing  or  originating  in  Tibet.  Finally,  Ni  and  Barazangi  (1983)  report  the 
results  of  a  detailed  study  in  India,  Tibet,  and  the  Himalayan  Arc  which  agrees  with 
other  results.  Chinn  et  al.  (1980)  studied  L,  and  SB  in  western  South  America  and 
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observed  that  Ls  propagation  was  efficient  only  parallel  to  the  strike  of  the  Andes. 
They  also  observed  conversion  of  Oceanic  SB  to  L(  in  areas  of  crustal  thickening. 
Isacks  and  Stephens  (1975)  also  observed  conversion  of  SB  to  L,  both  in  earthquakes 
from  the  West  Indies  observed  in  Eastern  North  America  and  in  earthquakes  near 
southernmost  Mexico  observed  on  the  northern  coast  of  the  Gulf  of  Mexico.  Some  of 
their  claims  of  observed  conversions  are  disputed  by  Shurbet  (1962,  1974,  1976).  Gre- 
gersen  (1978)  discusses  conversion  between  different  modes  of  Love  waves  and  between 
Love  and  Rayleigh  waves  at  an  ocean  continent  boundary.  These  studies  illustrate 


the  importance  of  a  clear  understanding  of  the  mechanisms  of  attenuation  of  L(  along 
mixed  paths.  Simply  stating  that  a  path  which  passes  under  an  ocean  but  is  not 
attenuated  implies  that  the  floor  of  that  ocean  has  a  continental  structure  could  be 
misleading.  A  clearer  understanding  of  the  effects  of  structural  transitions  on  the 
attenuation  of  Lf  is  important. 

Another  major  use  of  L(  waves  is  the  determination  of  magnitudes  of  explosions 
and  earthquakes  and  of  regional  attenuation  properties.  These  magnitudes,  usually 
referred  to  as  mbLjJ  and  other  types  of  magnitudes  are  compared  to  discriminate 
between  the  two  types  of  sources  (Blandford,  1982,  and  Pomeroy  et  al.  1982).  The 
ratio  of  body  wave  magnitude,  m*,  to  surface  wave  magnitude,  Ms,  is  one  such 
discriminant  whose  physical  basis  is  explained  by  Stephens  and  Day  (1985).  In  order 
to  examine  whether  this  type  of  discriminant  can  be  generalized  to  periods  shorter 
than  than  those  considered  when  using  the  twenty  second  surface  wave  magnitude, 
new  empirical  relations  for  short  period  phases  such  as  Lt  (Baker,  1970)  and  R(  (Bas¬ 
ham,  1971)  were  derived  to  measure  short  period  surface  wave  magnitudes.  Nuttli 
(1973)  used  observed  attenuation  of  short  period  Rayleigh  waves  (L(,  and  Rg)  to 
present  a  theoretical  explanation  and  generalization  of  the  empirical  relations  and 
defined  their  limitations.  He  explained  how  to  use  data  of  L,  amplitude  versus 


distance  to  determine  the  value  of  7,  the  coefficient  of  anelastic  attenuation.  Nuttli 
(1978)  extended  the  definitions  from  one  second  waves  to  frequencies  as  high  as  ten 
Hertz.  These  studies  provided  a  more  quantitative  method  for  mapping  the  attenua¬ 
tion  of  Lg  than  simply  recording  whether  Lg  was  observed,  not  observed,  or 
attenuated.  Thus,  a  series  of  studies  mapping  the  attenuation  of  Lg  using  the  determi¬ 
nation  of  7  were  initiated.  An  extensive  study  of  the  United  States  at  0.5  to  2  Hz.  had 
previously  been  completed  by  Sutton  et  al.  (1967)  using  other  methods.  Nuttli  (1973) 
developed  an  Lg  magnitude  scale  for  the  central  United  States.  Herrmann  and  Nuttli 
(1982)  studied  the  relation  of  this  scale,  m^  to  the  local  magnitude  scale,  ML.  The 
attenuation  observed  was  used  to  model  ground  motions  at  regional  distances  due  to 
earthquakes  and  assess  the  possible  destructiveness  of  the  L(  part  of  the  wavetrain 
(Herrmann  and  Nuttli,  1975a,  1975b).  An  additional  study  of  this  type  of  problem 
using  added  data  was  done  by  Street  (1984).  Street  et  aL  (1975)  studied  Lg  spectra  of 
United  States  earthquakes.  The  results  were  the  basis  of  a  study  of  Lg  moments,  mag¬ 
nitudes,  and  intensities  (Street  and  Turcotte,  1977)  and  the  development  of  an  Lg 
magnitude  scale  analogous  to  that  of  Nuttli  (1973)  for  the  northeastern  United  States 
and  southeastern  Canada  (Street,  1976).  Jones  et  al.  (1977)  studied  attenuation  of  Lg 
in  the  southeastern  United  States.  Nuttli  (1980)  considered  a  region  including  Iran. 
Bollinger  (1979)  again  studied  the  northeastern  United  States  obtaining  results  con¬ 
sistent  with  Street  (1976).  Nuttli  (1981)  examined  attenuation  of  Lg  in  western  and 
central  Asia  for  both  earthquakes  and  explosions.  Barker  et  al.  (1981)  studied  the 
crustal  structure  at  the  Nevada  Test  Site  using  Lg  and  Pg.  Chung  and  Bernreuter 
(1981)  reviewed  the  studies  of  regional  relations  between  various  magnitude  scales 
including  m^(.  Nicolas  et  al.  (1982)  studied  attenuation  of  regional  phases  in  western 
Europe.  Gupta  et  al.  (1980)  considered  both  earthquakes  and  explosions  in  western 
Russia,  and  Gupta  et  al.  (1982)  studied  variations  in  horizontal  to  vertical  Lg 
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amplitude  ratios.  Ebel  (1982)  criticizes  the  m^  scale  as  being  inappropriate  at  five  to 
ten  Hertz  for  earthquakes  in  the  northeastern  United  States.  Dwyer  et  al.  (1983) 
presented  a  frequency  dependent  study  of  L(  in  the  central  United  States  which 
addresses  this  criticism  of  the  m^  scale.  Singh  and  Herrmann  (1983)  present  a 
regionalized  map  of  Lg  attenuation  in  the  continental  United  States  for  periods 
between  0.5  and  3.5  Hz.  Herrmann  and  Kijko  (1983a,  1983b)  also  investigate  neces¬ 
sary  changes  to  the  L,  magnitude  scale  to  allow  use  of  data  for  instruments  with 
predominant  periods  other  than  the  one  second  period  of  the  WWSSN  instrument. 
Campillo  et  al.  (1984)  studied  the  excitation  of  spectra  and  geometrical  attenuation  of 
regional  phases  including  Lg.  Der  et  al.  (1984)  studied  the  coherence  of  the  L(  wave 
crossing  a  given  region.  All  these  studies  are  based  on  the  interpretation  of  Lg  as  a 
superposition  of  higher  mode  Love  and  Rayleigh  waves.  They  all  assume  that  the  L( 
waves  are  traveling  through  a  relatively  constant  thickness  crustal  layer,  and  most  do 
not  account  for  reflection  of  diffraction  of  L,  energy  from  changes  in  thickness  of  the 
crust.  Clearly,  it  is  important  to  understand  if  such  changes  will  produce  significant 
effects  not  accounted  for  in  the  interpretations  given. 

Many  attempts  to  understand  the  propagation  of  seismic  disturbances  across 
regions  of  varying  structure  such  as  transition  zones  have  been  made.  First,  simple 
models  were  used  and  analytic  solutions  were  derived  for  soluble  special  cases,  then 
increasingly  complicated  models  were  considered  as  available  computational  power 
increased.  The  types  of  models  that  have  been  used  to  approximate  transition  regions 
can  be  separated  into  several  types  which  are  illustrated  in  Figure  1. 

The  simplest  type  of  model  (Figure  la)  consists  of  two  homogeneous  layered 
regions,  11’,  and  3  3’,  with  layer  over  half-space  structures  separated  by  a  vertical 
boundary  or  a  region  2  2’  in  which  density,  p,  rigidity  p,  and  SH  wave  velocity  ,  0, 
vary  smoothly  between  their  values  in  the  regions  1  1’  and  3  3’.  Sato  (1961a)  discusses 
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Figure  1:  Types  of  models  used  in  studies  of  propagation  of  Love  waves  across 
continent  ocean  boundaries.  Part  a)  shows  the  simplest  model,  two  regions  labeled  1 
and  3  separated  by  a  vertical  boundary  or  a  intermediate  region,  2,  of  varying  proper¬ 
ties.  Part  b)  shows  another  type  of  simple  model,  a  layer  of  a  half-space  with  a  sudden 
step  change  in  the  thickness  of  the  layer.  Part  c)  shows  a  more  genera)  type  of 
simplified  model  with  a  more  gradual  change  in  thickness  rather  than  a  step.  Part  d) 
shows  the  type  of  model  used  in  this  study,  which  includes  smooth  thickness  variation 
both  at  the  surface  and  the  base  of  the  crustal  layer. 
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to  assure  the  transmission  of  fundamental  mode  Love  waves  across  the  boundary 
without  modal  conversion.  He  derives  analytic  expressions  for  the  transmitted  and 
reflected  waves,  the  phase  and  group  velocities,  and  the  reflection  and  transmission 
coefficients  in  each  case.  The  case  for  which  no  transition  regions  is  present  has  since 
been  used  by  Alsop  (1966),  Boore  (1970),  and  Gregersen  and  Alsop  (1974,  1976)  as  a 
test  case  for  more  generalizable  numerical  methods.  More  recently,  Kennett  (1973)  has 
developed  a  numerical  technique  to  consider  the  problem  of  seismic  waves  interacting 
with  a  layer  or  layers  in  which  properties  change  across  a  surface  perpendicular  to  or 
at  a  specified  angle  from  the  layering. 

The  next  level  of  complexity  that  can  be  introduced  into  the  model  is  the  inclu¬ 
sion  of  a  change  in  the  thickness  of  the  surface  layer  at  the  discontinuity  between  the 
two  structures.  This  type  of  model  is  illustrated  in  Figure  lb.  In  most  cases  the 
material  properties  do  not  change  when  the  layer  thickness  changes.  Sato  (1961b) 
used  the  Wiener-Hopft  technique  to  obtain  analytic  expressions  for  the  transmitted 
and  reflected  waves  for  a  model  with  a  surface  step  with  height  h  much  less  than  the 
wavelength  of  the  incident  wavefield  X.  Using  these  solutions  he  obtained  approxi¬ 
mate  expressions  for  energy  reflection  and  transmission  coefficients  for  a  surface  step 
model  using  a  Green’s  function  method.  Hudson  and  Knopoff  (1964)  calculated 
expression  for  the  motion  and  the  reflection  and  transmission  coefficients  for  the  sur¬ 
face  step  model.  They  did  not  need  or  use  the  h<<\  approximation.  Alsop  (1966) 
developed  an  approximate  method  for  determining  Love  wave  transmission  and 
reflection  coefficients  in  a  model  of  type  lb.  This  method  assumes  that  all  energy 
remains  in  Love  waves,  introducing  errors,  particularly  at  intermediate  periods,  if  any 


energy  is  converted  between  modes  or  into  forms  other  than  modal  energy.  The  latter 
two  methods  had  the  advantage  that  higher  modes  could  be  used  and  that  no  specific 
interrelation  of  elastic  properties  between  the  different  layers  or  regions  was  necessary. 
Gregersen  and  Alsop  (1974,  1976)  extended  the  method  of  Alsop  (1966,1968)  to  the 
determination  of  transmission  coefficients  for  the  case  of  non-normally  incident  Love 
waves.  At  normal  incidence  results  correspond  well  to  those  found  in  previous  studies, 
for  oblique  incidence  at  angles  less  than  forty  degrees  normal  incidence  is  a  good 
approximation.  Bose  (1975)  solves  the  surface  step  problem  using  an  integral  equation 
formulation  which  yields  asymptotic  solutions  for  large  distances  from  the  step.  His 
results  show  increases  in  amplitude  for  a  step  down  and  decrease  in  amplitude  for  a 
step  up,  consistent  with  results  derived  for  the  step  models  in  this  study.  Kazi  (1978a, 
1978b)  uses  the  variational  method  of  Schwinger-Levine  to  derive  solutions  to  the  sur¬ 
face  step  problem  that  account  for  the  Love  waves  converted  to  body  waves  at  the 
surface  step.  He  determines  the  reflection  and  transmission  coefficients  by  way  of  a 
scattering  matrix  and  evaluates  them  for  the  same  cases  as  Alsop  (1966)  and  Knopoff 
and  Hudson  (1964)  showing  that  the  transmission  coefficients  increased  after  the  cutoff 
frequency  rather  than  decreasing  as  was  observed  in  earler  studies  that  neglected  the 
body  waves.  Lapwood,  Hudson,  and  Kembhavi  (1973,  1975)  and  Lapwood  and  Hud¬ 
son  (1975)  used  a  similar  variational  technique  to  study  normally  and  obliquely 
incident  plane  waves  incident  on  a  layer  between  two  uniform  half-spaces.  Hudson 
(1977)  extended  their  method  to  the  surface  step  problem  but  concluded  that  the 
numerical  implementation  of  the  method  would  not  be  more  efficient  than  existing 
techniques.  Martel  (1980)  used  a  finite  element  technique  to  evaluate  propagation  of 
Love  waves  across  a  Moho  step.  Spatial  filtering  of  the  transmitted  and  reflected 
modes  to  determine  transmission  and  reflection  coefficients  allowed  the  isolation  of  the 
diffracted  body  wave  component. 


Studies  using  the  surface  step  model  and  the  Moho  step  model  have  been  general¬ 
ized  in  several  ways.  First,  and  most  relevant  to  this  study,  the  step  can  be  replaced 
by  a  transition  region  as  discussed  in  the  next  paragraph.  Other  types  of  generaliza¬ 
tion  allow  the  study  of  the  effect  of  topography,  of  a  short  section  if  thinner  crust 
referred  to  as  a  crustal  ’pinch’,  or  of  continuous  but  small  variations  in  crustal  thick¬ 
ness  on  the  incident  wavefield.  Gilbert  and  Knopoff  (1960)  used  a  method  similar  to  a 
Born  approximation  to  study  seismic  scattering  from  topography  and  to  produce 
approximations  valid  for  pulses  with  duration  short  with  respect  to  the  horizontal 
travel  time  across  the  irregularity.  Hudson  (1967)  extended  the  method  of  Gilbert  and 
Knopoff  (1960)  to  variations  in  elastic  properties  at  the  surface  such  as  sediment  filled 
depressions,  as  well  as  surface  distortions.  The  method  can  be  used  to  obtain  a  rough 
estimate  of  the  dimension  of  the  scatterer.  Boore  et  al.  (1971)  and  Boore  (1972)  study 
the  effect  of  topography  on  SH  waves  using  the  Finite  Difference  method  to  propagate 
a  simple  pube  through  the  structure.  This  approach,  unlike  the  previously  discussed 
methods,  allows  the  consideration  of  steep  topography  or  topography  with  wavelength 
similar  to  or  smaller  than  that  of  the  incident  waves.  Herrera  (1964a,  1964b) 
developed  a  perturbation  method  to  study  the  effects  on  a  propagating  seismic  distur¬ 
bance  of  a  crustal  layer  with  slowly  varying  thickness.  Mai  and  Herrera  (1965)  studied 
the  effect  of  a  short  section  of  thinner  crust,  where  the  difference  in  crustal  thickness 
between  the  thin  and  the  thick  crust  is  small  compared  with  the  wavelength  of  the 
incident  energy. 

The  next  increase  in  the  complexity  of  the  transition  model  is  the  introduction  of 
a  transition  region  in  which  the  crustal  thickness  varies  smoothly  between  the  crustal 
thicknesses  at  its  ends.  The  simplest  modeb  of  thb  type  are  shown  in  Sgure  1c.  In 
these  modeb  either  the  surface  or  the  Moho  have  a  smooth  slope  in  the  transition 
region.  Knopoff  and  Mai  (1967),  and  Knopoff  et  al.  (1970)  explained  an'  analytic 
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solution  for  this  type  of  model  when  the  slope  of  the  surface  (or  Moho)  of  the  transi¬ 
tion  region  is  small.  They  also  assumed  that  the  local  phase  velocity  was  the  same  as 
that  in  a  layered  structure  with  the  local  layer  thickness.  They  found  forward  scatter¬ 
ing  to  be  much  less  important  than  back  scattering.  Pec  (1967)  calculated  the  disper¬ 
sion  of  Love  waves  propagating  in  a  wedge  shaped  layer  such  as  the  transition  region 
in  models  of  type  lc.  He  considered  Love  waves  in  terms  of  the  constructive  interfer¬ 
ence  of  multiply  reflected  SH  rays  and  found  that  the  inclination  of  the  Moho  caused 
the  largest  changes  in  phase  velocity  and  amplitudes  at  short  periods.  Boore(1970) 
studied  the  propagation  of  a  simple  low  frequency  Love  wave  across  the  type  of  struc¬ 
ture  shown  in  figure  lc.  He  specified  the  displacements  at  t=0  and  t— dt  throughout 
the  grid  such  that  the  displacement  at  t=0  at  the  surface  of  the  layer  over  a  half¬ 
space  portion  of  the  model  was  a  Ricker  wavelet.  The  Finite  Difference  method  was 
then  used  to  calculate  phase  velocities  and  transmission  coefficients  at  various  points 
along  the  gTid  surface  for  a  one  hundred  twenty  kilometer  long  transition  region. 
Periods  of  twenty  to  one  hundred  thirty  seconds  were  used.  Phase  velocities  were 
lower  for  propagation  from  continent  to  ocean  than  from  ocean  to  continent.  He 
noted  that  in  the  region  of  the  transition  that  mode  conversions  and  conversion  to 
other  types  of  waves  seemed  to  be  important.  Lysmer  and  Drake  (1971,  1972)  and 
Drake(1972),  who  discusses  mainly  Rayleigh  wave  results  but  states  that  they  also 
apply  to  Love  waves,  use  a  Finite  Element  method  based  on  Zienkiewicz  and  Cheung 
(1967)  which  includes  a  rigid  grid  bottom  and  thus  allows  no  energy  to  escape  the 
grid.  The  formulation  also  requires  that  the  incident  modal  energy  is  exactly  equal  to 
the  sum  of  the  reflected  and  transmitted  modal  energy,  thus  disallowing  conversion  to 
other  types  of  waves.  Lysmer  and  Drake  (1971)  use  this  method  to  study  the  effect  of 
a  transition  of  type  lc,  with  length  one  hundred  twenty  kilometers,  or  Id,  with  a 
length  of  one  hundred  kilometers,  on  the  incident  fundamental  mode  Love  wave 
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energy.  With  their  constraint  no  difference  in  continent  to  ocean  and  ocean  to  con¬ 
tinent  phase  velocities  is  seen,  so  they  attribute  the  differences  to  body  wave  interfer¬ 
ence.  They  estimate  that  about  one  percent  of  the  energy  at  25s  period  is  transmitted 
without  mode  conversion  and  about  35%  is  transmitted  in  other  modes.  Drake(1971) 
discusses  mainly  Rayleigh  wave  results  but  states  that  they  also  apply  to  Love  waves. 
Drake  and  Bolt  (1980)  used  the  method  of  Lysmer  and  Drake  (1972)  to  study  a  more 
complicated  model  with  the  upper  and  lower  transition  surfaces  having  different 
lengths  to  model  phase  velocity  data  of  events  normally  incident  on  the  California 
continental  margin.  They  considered  periods  between  4.4  and  60  seconds  and  modeled 
only  incident  fundamental  mode  waves.  They  discuss  the  period  ranges  over  which 
fundamental  mode  transmission  is  small  and  conversion  to  particular  higher  modes  is 
large.  They  conclude  that  the  ocean  continent  boundary  strongly  increases  the 
attenuation  of  fundamental  mode  Love  waves.  Schlue  (1979)  shows  some  simple 
examples  of  a  three  dimensional  Finite  Element  solution  for  Love  wave  propagation. 

All  the  studies  discussed  in  the  previous  paragraph  used  periods  much  longer 
than  those  that  will  be  considered  in  the  following  discussions.  The  shorter  periods 
used  in  this  study  allow  the  examination  of  the  effects  of  transition  regions  with  length 
many  times  the  wavelength  of  the  incident  energy.  Most  previous  studies  considered 
transition  regions  with  lengths  comparable  the  wavelength  of  the  incident  energy.  The 
method  of  Lysmer  and  Drake  (1972)  purposely  excluded  the  possibility  of  conversion 
from  modal  energy  to  body  waves,  and  thus  excludes  the  possible  escape  of  energy 
from  the  system  as  these  body  waves  travel  out  the  bottom  of  the  grid.  This  escaping 
energy  is  shown  in  this  study  to  be  an  important  component  of  the  explanation  of  the 
attenuation  of  the  L(  phase  traveling  on  partially  oceanic  paths.  The  studies  using  the 
method  of  Lysmer  and  Drake  (1972)  used  periods  shorter  than  other  studies  but  longer 
than  those  used  here.  They  also  considered  only  fundamental  mode  Love  wave  input 


and  calculated  each  mode  separately,  while  the  driving  functions  used  in  this  study  are 
a  sum  over  a  range  of  frequencies  on  the  fundamental  and  first  five  higher  branches 
In  fact  most  of  the  earlier  studies  with  Finite  Element  or  Finite  Difference  considered 
the  fundamental  mode  at  one  frequency  rather  than  over  a  range  of  frequencies.  None 
of  these  studies  has  generated  realistic  seismograms,  but  instead  they  all  concentrated 
on  measuring  phase  velocities  and  transmission  and  reflection  coefficients.  In  light  of 
these  facts  it  is  clear  that  the  results  discussed  in  this  thesis  can  give  a  significant 
addition  to  the  understanding  of  the  propagation  of  Lf  waves  across  transition  regions. 

Finally,  the  content  of  each  chapter  will  be  summarized.  The  first  four  sections  of 
chapter  one  present  the  basics  of  and  the  notation  associated  with  the  Propagator 
Matrix  and  FE  techniques.'  Sections  five  and  six  discuss  the  theory  of  the  two  cou¬ 
pling  methods.  The  remainder  of  the  chapter  is  a  series  of  tests  of  the  accuracy  of  the 
coupling  technique  discussed  in  section  five  and  used  in  chapter  two.  Sections  one  to 
three  of  chapter  two  discuss  the  design  of  the  numerical  experiments  presented  in  the 
remainder  of  that  chapter,  and  in  section  two  of  chapter  three.  The  remainder  of 
chapter  two  discussed  the  effects  of  transition  of  L,  wavetrains  and  the  dependence  of 
those  effects  on  transition  length.  Chapter  three  first  presents  FE  results  investigating 
the  effect  of  the  length  of  the  oceanic  path  on  the  L,  wavefieid.  Then  examples  of  the 
accuracy  and  efficiency  of  the  coupling  method  discussed  in  section  six  of  chapter  one 
are  presented.  Finally,  this  method  is  used  to  determine  seismograms  after  propaga¬ 
tion  through  a  transition  and  then  oceanic  paths  of  different  lengths.  Reasons  for  the 
attenuation  of  L.  for  oceanic  path  lengths  longer  than  100-200  kilometers  are  dis- 
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Chapter  1 

Representation  Theorem  Coupling  of  the  Finite  Element 
and  Modal  Propagator  Matrix  Methods 


Intoduetion 

In  this  chapter  the  methods  for  Representation  Theorem  (RT)  coupling  of  finite 
element  (FE)  or  finite  difference  calculations  and  Harkrider’s  (Harkrider  1904,  1970) 
propagator  matrix  method  calculations  to  produce  a  hybrid  method  for  propagation  of 
SH  mode  sum  seismograms  across  paths  that  contain  regions  of  non  plane-layered 
structures  are  explained  and  developed.  The  FE  method  used  in  this  study  is  an 
extensively  modified  version  of  the  Stress  Waves  in  Solids  code  (Frazier,  Alexander  and 
Petersen,  1973).  The  coupling  methods  explained  in  detail  in  this  chapter  use  a  2-D 
Cartesian  FE  formulation.  Analogous  methods  for  the  3-0  method  follow  directly. 
The  hybrid  method  linking  propagator  matrix  results  to  FE  calculations  assumes  that 
the  wavefield  arrives  at  the  edge  of  the  FE  grid  after  propagating  many  wavelengths 
in  the  layered  structure.  Extensive  tests  illustrating  the  validity  and  accuracy  of  the 
implementation  of  this  method  are  presented  in  this  chapter.  These  tests  use  a  struc¬ 
ture  consisting  of  a  layer  over  a  half-space.  This  simple  structure  allows  the  calcula¬ 
tion  of  synthetic  seismograms  using  only  the  propagator  matrix  technique.  These  syn¬ 
thetic  seismograms  can  be  directly  compared  to  the  hybrid  results  to  ascertain  the 
accuracy  and  the  limitations  of  the  coupling  technique.  After  propagation  through  a 
path  of  length  much  less  than  the  source  to  FE  gTid  distance  in  the  FE  grid  the 
wavefields  can  be  coupled  back  into  a  layered  medium.  The  layering  in  the  second 
layered  medium  need  not  be  identical  to  that  in  the  first  layered  medium.  The  tests 
of  coupling  from  FE  to  propagator  matrix  are  discussed  in  chapter  3.  This  sequence  of 


procedures  can  be  repeated  any  number  of  times,  thus  including  several  regions  of 
complexity  in  the  source  to  receiver  path. 

These  hybrid  techniques  are  developed  to  study  the  propagation  of  surface  waves 
across  regional  transition  zones  or  other  heterogeneities  that  exist  in  part  of  a  longer, 
mostly  plane-layered,  path.  Examples  of  structures  of  interest  through  which  surface 
waves  can  be  propagated  using  these  techniques  include,  regions  of  crustal  thickening 
or  thinning  such  as  continent  ocean  transitions  or  basins,  anomalous  bodies  of  any 
shape  located  in  the  path,  topography,  and  sudden  transitions  from  one  structure  to 
another.  In  fact,  any  arbitrary  structure  of  interest  can  be  placed  in  the  FE  portion  of 
the  path. 

To  clearly  explain  the  methods  used  to  couple  finite  element  (FE)  and  propagator 
matrix  methods,  enabling  transmission  of  a  disturbance  along  part  of  its  path  from  the 
source  to  the  receiver  by  each  method,  it  is  useful  to  first  discuss  some  of  the  founda¬ 
tions  on  which  each  method  is  based,  and  some  of  the  details  of  the  implementation  of 
each  method.  The  propagator  matrix  technique  is  used  to  transmit  mode  sum  syn¬ 
thetics  through  sections  of  their  paths  that  consist  of  simple  plane  layered  structures. 
The  theory  used  to  derive  this  technique  is  explained  in  detail  by  Harkrider  (Harkrider 
1964).  His  numerical  implementation  of  these  ideas  forms  the  basis  for  the  slightly 
modified  codes  used  in  this  study.  The  fundamentals  of  this  method  which  are  central 
to  the  understanding  of  the  Representation  Theorem  coupling  techniques  and  the  gen¬ 
eration  of  the  mode  sum  synthetics  and  the  Green’s  Function  synthetics  used  in  these 
techniques  will  be  discussed  below.  The  modifications  to  the  numerical  implementa¬ 
tion  of  Harkrider  will  also  be  discussed.  The  FE  calculations  discussed  here  are  done 
using  an  extensively  modified  version  of  the  Stress  Waves  in  Solids  (SWIS)  code 
developed  by  Frazier,  Alexander,  and  Petersen  (Frazier  et  al.  1973,  Frazier  and  Peter¬ 
sen  1974).  The  basic  structure  of  SWTS  code  is  retained,  but  many  options  in  the  code 


are  removed,  and  the  detailed  implementation  is  changed  in  order  to  produce  a  large 
increase  in  the  speed  of  execution.  Further  options  not  available  in  the  version  of  the 
SWIS  code  used  were  added  to  enable  grid  points  to  be  constrained  to  move  with  a 
given  time  history,  to  allow  the  production  of  time  slices,  and  to  implement  hourglass 
corrections  for  the  SH  case.  Details  of  these  changes  and  their  impact  will  be  dis¬ 
cussed  later,  as  will  the  basic  formulation  of  the  SWIS  FE  method. 

The  Propagator  Matrix  Technique 

The  basic  idea  behind  the  propagator  matrix  method  as  implemented  by  Har- 
krider  is  to  transmit  the  disturbance  produced  by  a  source  within  a  layered  half-space 
structure  through  that  structure  by  combining  terms  that  describe  the  source,  the 
medium  response,  and  the  propagation  effects.  The  terms  that  apply  the  effects  for 
the  propagation  path  in  the  z  direction  are  in  the  form  of  propagator  matrices.  For 
source  and  receiver  both  at  depth,  they  are  separated  into  two  parts,  the  propagation 
in  z  from  the  source  to  the  surface,  and  the  propagation  in  z  from  the  surface  to  the 
receiver.  An  additional  propagation  term  adds  the  effect  of  the  propagation  in  the  r 
coordinate.  In  most  cases  the  general  form  of  a  modal  displacement  at  a  receiver  at 
depth  is 
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where  S  is  a  function  of  the  source  strength  and  geometry,  At  represents  the  medium 
response  for  a  surface  source  and  a  surface  receiver,  P  expresses  the  propagation  effects 
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in  direction  r, 
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depth  to  the  surface  which  because  of  reciprocity  can  be  expressed  as  a  modal  propap 
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surface  to  the  receiver  depth.  The  subscript  H  on  the  terms  which  transmit  the  modal 
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disturbances  in  z  denotes  homogeneous,  that  is  independent  of  and  not  containing  a 


source.  The  expressions  within  the  square  brackets  of  these  terms  represent  matrix 


quantities,  not  simple  ratios.  If  the  source  is  a  stress  source  rather  than  a  displacement 


source  then 
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depth  of  the  source.  If  stress  rather  than  displacement  is  to  be  recorded  at  the 
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receiver,  then  -  is  replaced  by  — —  — .  Iq  this  expression  is  the 
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rigidity  at  the  depth  of  the  receiver.  As  an  example  of  a  specific  case  of  this  type  of 


representation,  the  displacement  expressions  for  a  double  couple  source  of  arbitrary 


orientation  will  be  stated  in  terms  of  the  propagator  formulation.  This  expression  will 


then  be  used  to  obtain  the  expressions  for  displacements  from  a  dip  slip  and  a  strike 


slip  source.  These  expressions  will  be  used  extensively  in  subsequent  discussions.  Fol¬ 


lowing  this  discussion  the  modifications  made  to  the  propagator  matrix  codes  will  be 


explained. 


The  displacement  for  an  arbitrary  double  couple  source  follows  directly  from  the 


expression  for  the  SH  displacement  at  the  free  surface  produced  by  a  double  couple 


source  of  arbitrary  orientation  at  depth  h  (  Harkrider  1964,  1970  ). 
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In  equations  (2),  (3),  and  (4),  X  is  the  strike  of  the  double  couple  source,  6  is  the  dip,  <t> 
is  the  azimuth  to  the  station,  0  is  the  SH  wave  velocity,  px  is  the  density  at  the  depth 
of  the  source,  h  is  the  source  depth,  w  is  the  frequency,  r  is  the  distance  from  the 
source  to  the  receiver,  cL  is  the  Love  wave  velocity,  &nd  «e  respectively  the  rigi¬ 
dity  at  the  source  depth  and  at  the  receiver  depth,  Al  is  the  medium  response  for  a 


surface  source  and  receiver, 
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are  the  terms  that  transmit 


the  source  disruption  in  z  from  the  source  to  the  surface  receiver,  and  the  term  con¬ 
taining  the  Hankel  function  is  the  propagation  term  in  r.  The  source  term  is  defined 
to  contain  all  the  angular  dependence  on  6,  X,  and  as  well  as  the  term 
and  the  factors  of  p.  The  second  equality  in  the  first  expression  of  equation  (3) 


assumes  a  step  moment,  that  is  R(w)  “  — 
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It  is  well  known  that  a  double  couple  of  arbitrary  orientation  can  be  expressed  in 
terms  of  a  linear  combination  of  double  couple  sources  of  three  types,  vertical  strike 
slip,  vertical  dip  slip  and  45°  dip  slip.  Thus,  results  for  these  three  fault  types  can  be 
added  to  produce  results  for  an  arbitrary  orientation,  removing  the  necessity  to  repeat 
the  entire  procedure  for  each  orientation  to  be  studied.  Substituting  (3)  and  (4)  into 
(2)  and  evaluating  the  expression  for  each  of  the  three  basic  fault  types  yields,  (5),  the 
expressions  used  to  determine  displacement  at  the  surface.  For  a  vertical  strike  slip 
fault  6  *  90°,  and  \=0#,  In  the  following  development  part  a  of  each  set  of  equa¬ 
tions  will  be  the  expression  for  the  strike  slip  fault.  For  a  vertical  dip  slip  fault 
6  —  90°,  and  X  =-  90°.  Part  b  of  each  set  of  equations  below  will  be  the  appropriate 
expression  for  the  vertical  dip  slip  fault.  For  a  45°  dip  slip  fault  6  =*  45°,  and 
X  am  90°.  Part  c  of  the  next  group  of  equations  is  for  a  45°  dip  slip  fault. 
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Examining  equations  (5)  shows  that,  for  SH  waves,  any  fault  geometry  can  be 
modeled  using  a  linear  combination  of  only  the  vertical  dip  slip  and  the  vertical  strike 
slip  faults.  This  follows  from  the  observation  that  the  vertical  dip  slip  and  the  45°  dip 
slip  expressions  differ  only  by  a  function  of  the  azimuth  For  the  remainder  of  this 
development  only  vertical  dip  slip  and  vertical  strike  slip  fault  types  will  be  con¬ 
sidered.  This  is  equivalent  to  setting  sin^  =  cos2^  ®  1,  and  sin2^  —  0,  that  is 
0  =  90°.  For  a  given  value  of  4  any  fault  geometry  can  be  reconstructed  by  reintro¬ 
ducing  the  appropriate  values  of  sin^  ,  cos2^  ,  and  sin24  into  the  coefficients  of  the 
linear  combination  of  the  remaining  two  factors.  The  displacements  determined  by 
equation  (5)  are  those  produced  at  the  surface  by  a  double  couple  source  at  depth  h. 
To  determine  the  displacements  at  depth,  a  term  propagating  the  energy  from  the 
surface  to  the  depth  of  the  receiver  must  be  included.  For  a  receiver  at  depth  z  the 
modal  displacement  expressions  become 
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Most  of  the  changes  made  to  the  propagator  matrix  codes  of  Harkrider  were 
minor.  The  only  major  changes  were  the  addition  of  an  option  to  evaluate  the 


analytic  expressions  for  the  stresses  and  the  implementation  of  an  alternate  code 
which  uses  the  line  source  rather  than  the  point  force,  double  couple,  or  explosive 
sources  already  available.  The  implementation  of  the  stress  calculation  will  be  dis* 
cussed  later.  The  construction  of  the  line  source  code  required  changes  to  the  source 
terms  and  the  modal  propagation  terms  that  propagate  in  r.  The  basic  form  of  the  z 
dependent  modal  propagator  terms  is  identical  regardless  of  the  source  used,  so  no 
changes  were  necessary  in  the  calculation  of 
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calculate  displacement,  stress,  and  displacement  and  stress  Green’s  functions.  The 
minor  changes  made  to  the  codes  include  the  restructuring  of  the  order  of  operations 
to  allow  the  efficient  calculation  of  particular  types  of  sets  of  seismograms.  To 
efficiently  calculate  sets  of  displacement  and  stress  seismograms  in  which  successive 
seismograms  are  at  the  same  horizontal  distance  r,  and  at  a  series  of  receiver  depths  at 
intervals  of  Az,  the  operations  were  reordered  to  avoid  unnecessary  recalculation  of 
intermediate  results.  The  propagation  terms  in  r  need  to  be  calculated  only  once,  as 
does  the  source  to  surface  term.  To  efficiently  calculate  sets  of  displacement  and 
stress  line  source  Green’s  functions  the  calculations  were  reordered  in  a  different 
manner  so  that  the  surface  to  receiver  propagator  term  and  the  propagation  term  in  r 
needed  to  be  evaluated  only  once.  The  codes  still  retain  the  ability  to  calculate 
seismograms  for  sets  of  receivers  or  Green’s  functions  for  sets  of  sources  which  do  not 
form  depth  sections.  However,  for  these  sets  the  order  of  operations  is  not  necessarily 
optimal  and  extra  time  may  be  taken  in  calculation  of  displacements  or  stresses. 

The  Modified  SWIS  Finite  Element  Code 


The  SWIS  FE  program  is  a  flexible  code  capable  of  calculating  the  propagation  of 
stress  waves  through  solids  of  arbitrary  complexity.  Properties  such  as  density,  and 
elastic  wave  velocity  are  defined  at  a  grid  of  nodes  whose  locations  can  be  defined  in 


curvilinear  coordinates.  These  properties  are  then  used  to  determine  the  location, 
stress,  and  motion  of  the  nodes  in  the  grid  as  a  function  of  time.  Next,  it  is  useful  to 
give  a  simple  step  by  step  outline  of  the  procedures  used  in  the  SWIS  FE  code. 

First,  several  quantities  are  defined  at  each  node  or  at  each  element  in  the  grid  to 
provide  the  initial  conditions.  In  particular,  the  properties  that  are  initialized  at  each 
node  before  the  first  time  step  are,  the  nodal  position  X(t),  nodal  displacement  li.(t), 

and  nodal  velocity  U  |t  -  j .  The  properties  defined  for  each  element  are  the  force 

as  a  function  of  time,  £(t),  resulting  from  body  forces  JJt)  within  the  element  and  from 
tractions  i(t)  applied  to  the  surface  of  the  element,  and  the  elastic  constants  as  deter¬ 
mined  by  specification  of  density  and  elastic  wave  velocities.  The  values,  within  a  ele¬ 
ment,  of  the  variables  initialized  on  each  node  are  defined  by  interpolation  between 
the  nodal  values  at  the  nodes  surrounding  that  element.  In  order  to  define  the  form  of 
the  interpolation  functions  the  coordinate  systems  which  define  the  location  of  a  point 
within  an  element  and  the  location  of  an  element  within  the  body  must  be  defined  and 
related.  The  ’natural’  or  local  coordinate  system  defines  the  location  of  any  point  in 
the  element  with  respect  to  the  nodes  surrounding  that  element.  Each  node  is  defined 
to  be  a  unit  distance,  in  each  coordinate,  from  the  element  center.  Let  be  the  coor¬ 
dinates  of  a  point  in  element  e  in  this  ’natural’  coordinate  system,  and  let  Z  ™  be  the 
coordinates  of  node  m,  one  of  the  nodes  bordering  element  e,  in  the  ’natural’  coordi¬ 
nate  system.  This  means  that  Z,m  is  always  ±1.  A  Cartesian  coordinate  system  is 
defined  within  the  volume  V  containing  the  body  composed  of  the  elements.  In  this 
Cartesian  coordinate  system,  the  ith  component  of  the  position  of  node  m  bordering 
element  e  will  be  denoted  X”,  and  the  itk  coordinate  of  a  point  z  within  element  e  will 
be  denoted  Xj.  The  Cartesian  coordinates  Xj  and  Xjm  can  be  transformed  to  an  arbi¬ 
trary  curvilinear  coordinate  system  y-,  and  Yin  using  the  well  known  relations  for  the 
transformation  of  quantities  between  Cartesian  and  general  curvilinear  coordinates. 
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Since  all  calculations  in  this  study  are  performed  in  Cartesian  coordinates,  the 
transformation  to  general  curvilinear  coordinates  will  be  omitted  in  further  discussion 
of  the  SWIS  FC  method.  Now  it  is  possible  to  define  the  interpolation  function  used 
to  determine  positions,  dislocations,  velocities,  accelerations  and  stresses  within  an  ele- 
ment.  In  ’natural’  coordinates  the  interpolation  function  used  to  define  the  coordi¬ 
nates  of  a  point  z  within  element  e  has  the  form 

p"W  =  •jr.nu  +  «*•)  p> 

where 

Pm(g)  the  interpolation  function  at  point  z  for  position  of  node  m  bordering 

element  e 

D  number  of  space  dimensions 

Using  the  interpolation  function  it  is  possible  to  define  the  value  of  any  of  the  quanti¬ 
ties  initially  defined  at  each  node  at  any  point  within  an  element,  in  terms  of  the 
Cartesian  coordinates  of  the  nodes  surrounding  the  element.  The  location  of  the 
evaluation  point  within  the  element  is  expressed  in  the  ’natural’  coordinate  system, 
For  the  position,  displacement,  applied  force,  or  velocity  at  any  point  in  an  element 
surrounded  by  2°  nodes  the  relation  is  or  is  analogous  to  equation  (9). 

xj  -  £  P"U)XiB  (8) 

n«l 

For  the  stresses  the  relation  is 


The  SWIS  code  is  based  upon  a  one  point  quadrature  integration  rule,  so  the  point  i 
used  to  evaluate  these  expressions  is  the  centroid  of  the  element.  Given  the  initial 
values  the  SWIS  code  uses  one  of  several  series  of  calculations  to  advance  the  proper¬ 
ties  of  the  nodes  of  the  grid  one  time  step,  At.  Only  the  option  used  for  the 


calculations  in  this  study  is  discussed  below.  To  update  the  values  of  the  initialized 
quantities  by  At,  thus  yielding  the  configuration  of  the  grid  at  time  t+At  the  follow¬ 
ing  time  centered  explicit  FE  integration  procedure  is  used; 

(1)  Compute  strain  and  thus  stress  at  the  centroid  of  each  element  m 

=  2*{Ii^-ui(t)  +  li^Mu^t)  j  4-  x*iji£Muk(t)  (10) 


(2)  Compute  the  restoring  forces  for  each  element  m 


(id 


Then  use  the  values  of  restoring  force  R^™  at  each  element  m  surrounding  a  given  node 
to  determine  a  net  restoring  force  Rj  at  each  node. 

(3)  Integrate  the  equation  of  motion  to  give  the  values  of  location,  stress  and  motion 
at  each  node  at  time  t  +  At. 


,  .  Atl  (F,W +  *,(•)) 

U,(t  +  At)  -  U,(t)  +  AtjU  |t  -  £5.  +  At - ^ - 


(12) 


The  basic  series  of  calculations  described  above  form  the  core  of  the  SWIS  code  and 
were  retained  in  principle  in  the  modified  and  accelerated  version  of  SWIS  used  in  this 
study. 

In  the  discussion  of  the  basics  of  the  FE  method  used  no  mention  was  made  of 
the  exact  method  by  which  a  source  can  be  introduced  into  the  FE  calculations.  In 
the  basic  version  of  SWIS  a  point  source  or  a  source  of  finite  spatial  dimension  can  be 
defined  by  constraining  a  node  or  nodes  in  the  grid  to  have  a  specified  displacement  at 
time  zero.  Alternately,  given  body  forces,  Fj,  can  be  applied  at  the  centroids  of  one  or 
more  elements  at  time  zero  to  create  a  source.  Clearly,  a  more  general  source  would 


be  desirable,  and  in  fact  necessary  for  the  RT  coupling  of  Propagator  Matrix  or  other 
results  into  a  FE  calculation.  To  couple  energy  from  a  source  outside  the  FE  grid  into 
the  FE  calculation  it  is  necessary  to  specify  the  displacement  time  history  of  a  column 
of  nodes.  An  additional  option  to  allow  the  specification  of  displacement  as  a  function 
of  time  for  any  node  or  nodes  within  the  grid  was  added  to  the  basic  SWIS  code.  This 
option  allows  the  specification  of  a  time  history  constraint  at  up  to  two  hundred 
nodes.  However,  if  time  history  constraints  are  specified  at  more  than  one  node,  they 
must  have  a  common  duration  at  all  nodes  used.  Constraints  used  to  define  boundary 
conditions  are  defined  separately  and  persist  for  the  entire  calculation.  The  time  his¬ 
tory  constraints  can  have  any  duration.  If  the  duration  of  the  constraint  is  less  than 
the  duration  of  the  calculations  then  the  constrained  nodes  are  released  when  the  time 
t  within  the  calculation  exceeds  the  duration  of  the  constraint.  If  the  duration  of  the 
constraint  exceeds  the  duration  used  to  define  the  length  of  the  calculation  then  only 
the  first  portion  of  the  constraining  displacement  time  histories  are  used. 

A  second  modification  to  the  SWIS  FE  code  was  added  to  enable  the  production 
of  time  slices,  or  ’snapshots’  of  the  motions  at  each  node  within  the  grid  at  a  chosen 
time  t.  This  option  allows  the  code  to  write  a  file  which  contains  the  displacements  at 
every  node  in  the  grid,  or  at  every  node  in  a  decimated  grid  every  n  time  steps.  If  a 
decimated  grid  is  chosen  it  is  specified  by  an  integer  decimation  factor.  If  this  factor 
has  a  value  1  then  in  a  given  row  (H)  nodes  are  skipped  between  each  node  for  which 
a  value  is  recorded,  and  (1-1)  rows  of  nodes  are  skipped  between  each  row  where  dis¬ 
placements  are  recorded.  These  displacements  are  then  plotted  to  produce  a  series  of 
figures  showing  the  displacement  as  a  function  of  location  within  the  grid  at  a  selec¬ 
tion  of  times.  The  sequence  of  plots  illustrates  how  the  energy  propagates  through  the 
grid  as  a  function  of  time. 
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Tbe  final,  and  most  important  addition  to  the  SWIS  code  was  the  SH  hourglass 
correction  terms.  The  finite  element  (FE)  code  uses  a  one-point  quadrature  integra¬ 
tion  rule  so  stresses  are  evaluated  only  at  the  center  of  mass  of  each  element.  Conse¬ 
quently  any  stress  state  that  produces  non-constant  stress  within  an  element  is  not 
treated  completely  in  the  calculation.  Such  a  stress  state  can  be  decomposed  into 
the  sum  of  a  constant  stress  state,  which  the  one  point  integration  represents  exactly, 
and  a  non-constant  stress  state  with  zero  stress  at  the  centroid,  known  as  the  hour¬ 
glass  mode.  To  account  for  the  hourglass  modes  in  a  code  which  uses  one  point 
integration  an  additional  term  is  added  to  the  restoring  force.  The  form  of  this  term 
for  SH  or  antiplane  motion  is  explained  and  derived  below.  For  the  in  plane  or  P-SV 
problem  corresponding  hourglass  modes  occur  as  discussed  by  Kosloff  and  Frazier 
(1978).  In  all  cases  the  hourglass  term  in  the  restoring  force  is  calculated  using  similar 
procedures. 

Figure  la)  shows  a  2-D  rectilinear  element  of  the  type  used  in  the  SH  Finite  Ele¬ 
ment  calculations,  and  the  form  of  the  hourglass  deformation  present  for  SH  motion  of 
such  an  element.  This  deformation  lifts  one  pair  of  diagonally  opposite  element 
corners  and  depresses  the  other  pair  of  corners  by  the  same  amount,  producing  no 
stress  or  dispacement  at  the  centroid  of  the  element.  Since  the  Finite  Element  code 
depends  upon  evaluation  of  stress  only  at  the  centroid  of  each  element,  it  is  singular 
with  respect  to  such  deformation.  A  correction  to  the  restoring  forces  must  be  made 
to  offset  this  hourglass  motion.  This  correction  must  be  singular  with  respect  to  con¬ 
stant  stress  states  and  rigid  body  translations.  The  constant-stress  states  produce 
motions  as  shown  in  Figure  lb)  and  lc).  It  is  clear  that  no  sum  of  these  motions  and 
a  translation  in  y  will  produce  the  deformation  shown  in  la).  Thus,  the  hourglass 
deformation  is  indeed  singular  with  respect  to  translations  and  constant-stress  states. 


Figure  1:  Types  of  SH  motion  of  a  2-D  rectilinear  element,  a)  shows  the  hourgla» 
mode  for  SH  motion.  Arrows  at  the  element  corners  show  the  direction  of  displace* 
ment  in  y.  The  displacement  and  stress  at  the  centroid  (coordinate  origin)  of  th«  ele¬ 
ment  are  zero,  b)  and  c)  show  the  motions  accounted  for  in  the  uncorrected  stiffness 
matrix,  the  constant  stress  modes. 


For  a  body  undergoing  pure  antiplane  (SH)  deformation,  perpendicular  to  the  x-z 


plane,  the  following  conditions  are  satisfied: 


») 

*») 

<0 

d) 
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'«7 


du  j 
dx 
du- 
~di 


(13) 


These  relations  are  used  as  a  basis  to  derive  the  form  of  this  hourglass  correction  fac¬ 
tor  to  the  restoring  forces  for  the  case  of  SH  motion.  Using  equations  (7),  (8)  and  (9) 
and  the  equation  (14)  below,  from  Zienkiewicz  (1971),  the  hourglass  correction  to  the 
restoring  force  at  node  m,  Rjm,  can  be  calculated.  The  details  of  this  calculation  are 
given  in  the  following  paragraphs.  Equation  (14)  expresses  the  hourglass  restoring 
forces  on  each  node  in  terms  of  the  interpolation  functions,  the  tractions  on  the  sur¬ 
faces  of  the  surrounding  elements  and  the  body  forces  in  those  elements. 


Rim  *  £/hipm(i)dV  +  E/P^TjdS  (14) 

««v.  *»sr 

where 

hj  body  forces  in  volume  V,  i  component  i«*»l,  D 

S,  the  portion  of  the  element  surface  S  on  which  tractions  Tj  are 

specified 

£  sum  over  all  elements  bounding  the  node 

•m 

To  calculate  the  hourglass  restoring  forces  the  traction,  Tj,  on  the  element  sur¬ 
faces  are  needed.  The  quantities  used  as  tractions  for  the  purpose  of  calculating  hour¬ 
glass  restoring  forces  are  the  stresses  9^  and  9tJ.  Consider  the  2-D  rectilinear  element 
illustrated  in  Figure  2,  with  a  linear  varying  traction  9 „  applied  to  its  sides  at  z=±b 
and  9gj  applied  to  its  sides  at  x*±a.  Let  the  form  of  the  tractions  9 ^  be 


Figure  2:  Configuration  of  an  element  need  to  calculate  hourglass  restoring  forces 
for  2-D  SH  motion.  The  diameters  of  the  circles  on  the  element  edges  indicate  the 
amplitudes  of  the  applied  tractions  Ti  (Tf  -•  x~±a,  T|  ■■  o,lyi  z— ±b)  at 

selected  points.  Variation  between  illustrated  points  is  linear.  Circles  with  embedded 
crosses  indicate  motion  into  the  page.  Those  with  embedded  dots  indicate  motion  out 
of  the  page.  Z{  are  the  natural  coordinates  of  the  nodes.  The  lengths  a,  b  are  the 
dimensions  of  the  element  in  the  Cartesian  coordinate  system. 


From  this  relation  uy  can  be  determined  by  integrating  <rxy  with  respect  to  x. 

uy  =  J<t0  zdx.  =  (TqZX  +  f(z)  (16) 


However,  at  x=0,  z=»b,  uy=*0  so  f(z)*=0,  Now  uy  at  each  of  the  j  nodes  surrounding 
element  e  may  be  determined.  Expressing  the  uy  for  node  m,  m=l,..,j,  as  the  mth  ele¬ 
ment  of  a  vector  {Uy},  and  defining  the  unit  vector  (1}  to  carry  the  sign  information 

as 


1,  1,  1 


-} 


(17) 


enables  the  y  displacements  at  the  nodes  to  be  written 

h 

Simple  manipulation  of  (18)  leads  to  the  following  expression  for  <r0 

(,9) 

It  should  also  be  noted  that  equations  (16),  (13a),  and  (13d)  lead  directly  to  an  expres¬ 
sion  for  ax j 

a ^  (20) 

Now  we  have  the  expression  for  the  traction  T-t  ***  axv  <rx]  to  be  used  in  equation 
(9)  to  calculate  the  restoring  forces.  For  a  2-D  quadrilateral  element  subject  to  (13) 
the  body  forces  hj  are  zero  so  the  first  term  in  (9)  vanishes  and  the  interpolation  func¬ 


KR-  1) 
M 1. 1) 
M-1,-1) 
M  i,-i) 


gpab 


W 


(18) 


tions  become 
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p-to-^O  +  Z‘z”l  (■  +  *A”I 


Substituting  (21)  into  (14)  gives  the  following  expressions  for  Rm,  the  hourglass  restor¬ 


ing  forces. 


%  b  /  2  2 

r1  “  -/t  ( l-*i )  *0*  dx  +  /4  (1+**)  dz  =  <r0  *  b- 


R2  —  /i.  (l+zx)  <ToX  dx  +  /I  (14**2)  *0*  dz  —  -*o  ^ -  4- 
2  -a  *  *» 


*  ^  1  k 

R3  ™  '/t  (1“zi)  *0*  dx  -  /l  (l-za)  <roZ  dz  —  -<r0  (22) 

2  -b  2  6 


%  t 

R<  ~  /t  (1+zi)  *0*  dx  -  /i  (1-Zj)  <roZ  dz  —  <r0  (*  ±  b *? 
-* 2  -b 2  J 


Thus,  in  vector  form  the  restoring  forces-are 


ETa{'l 


Substituting  (15)  in  (19)  gives  the  value  of  the  restoring  forces  used  in  the  SH  FE 


code.  Square  elements  were  used,  that  is  a*b,  so  the  final  form  of  the  hourglass  res¬ 


toring  force  corrections  are: 


MM  M 


This  restoring  force  has  been  implemented  in  the  version  of  SWIS  used  for  all  subse¬ 


quent  calculations 


It  remains  to  be  shown  that  these  corrections  are  necessary  and  adequate  to  allow 


the  FE  code  to  produce  accurate  results.  To  illustrate  the  effect  of  the  corrections  two 


test  examples  were  completed,  each  of  them  with  and  without  the  hourglass  restoring 


-81- 


force  terms.  The  first  example  uses  a  point  source  located  at  one  of  the  nodes  within 
the  grid.  This  source  configuration  is  known  to  be  a  severe  test  of  the  stability  of  a 
FE  or  finite  difference  (FD)  calculation.  The  second  example  uses  a  two  step  pro¬ 
cedure,  the  validity  of  which  will  be  demonstrated  later,  to  couple  energy  propagated 
from  a  distant  line  source  into  the  FE  calculation  through  the  edge  of  the  FE  grid. 
This  procedure  produces  a  more  stable  source  and  introduces  hourglass  instabilities 
only  in  the  FE  portion  of  the  path.  Consequently,  smaller  differences  are  introduced 
into  the  results  when  hourglass  restoring  forces  are  omitted. 

Results  of  the  first  set  of  calculations,  using  the  line  source  within  the  grid,  are 
shown  in  Figure  3.  The  line  source  is  located  at  a  depth  of  10  km  from  the  free  sur¬ 
face  and  a  horizontal  distance  of  fifty  kilometers  from  the  left  hand  grid  edge.  The 
grid  has  dimensions  of  200x100  nodes,  and  node  spacing  of  dx»dy®=0.5  km.  All  cal¬ 
culations  use  a  time  step  of  0.05s.  The  displacement  time  history  of  the  source  is  a 
triangle  eight  time  steps  wide  with  a  rise  time  of  four  time  steps.  Each  group  of  three 
seismograms  shown  in  Figure  3  is  recorded  at  a  given  node.  The  nine  nodes  for  which 
seismograms  are  illustrated  surround  four  adjacent  elements  that  share  a  common 
node.  That  common  node  is  12.5  km  from  the  source,  at  the  same  depth  as  the 
source.  The  geometric  arrangement  of  the  groups  of  seismograms  in  Figure  3  is  identi¬ 
cal  to  that  of  the  nodes  at  which  they  are  recorded.  The  maximum  amplitude  of  each 
trace  is  normalized  to  one.  The  uppermost  seismograms  in  each  group  of  three  seismo¬ 
grams  in  Figure  3  is  calculated  by  direct  numerical  evaluation  of  the  analytic  expres¬ 
sion  for  the  propagation  of  waves  from  a  line  source  through  a  half-space.  This 
expression  is 


t 


t 


Figure  3:  Waveforms  for  synthetic  and  FE  seismograms  calculated  at  nine  adja¬ 
cent  nodes.  Amplitude  of  each  trace  is  normalized  to  one.  Each  group  of  three  traces 
is  recorded  at  one  node.  The  groups  of  traces  are  displayed  in  the  same  geometrical 
arrangement  as  the  nodes  at  which  they  are  recorded.  The  central  node  is  at  a  depth 
of  10  km,  12.5  km  from  a  10  km  deep  source.  The  spacing  between  the  illustrated 
nodes  is  0.5  km.  In  each  group  the  top  trace  is  synthetic,  the  second  trace  is  a  FE 
result  including  hourglass  restoring  forces,  and  the  third  trace  is  a  FE  result  without 
hourglass  restoring  forces. 


Where  rj  is  the  direct  source  to  receiver  distance,  r2  is  the  source  to  receiver  distance 
for  the  reflection  from  the  free  surface,  v  is  the  SH  wave  velocity  in  the  half-space,  t  is 
the  time  of  observation  with  t=-0  being  the  origin  time,  and  x  and  z  are  the  receiver 
coordinates.  The  first  term  on  the  right  hand  side  of  this  equation  gives  the  direct 
arrival,  the  second  term  the  reflection  from  the  free  surface.  The  first  peak  in  each 
synthetic  is  the  direct  arrival,  the  second  peak  is  the  free  surface  reflection.  The  cen¬ 
tral  seismogram  in  each  group  is  the  FE  result  including  the  hourglass  restoring  forces. 
The  last  seismogram  in  each  group  is  the  FE  solution  omitting  the  hourglass  terms. 
The  waveforms  and  amplitudes  for  the  first  two  traces  in  each  group  differ  mainly  in 
that  some  spurious  oscillations  are  present  in  the  FE  solution.  These  oscillations  are 
due  to  instabilities  caused  by  the  finiteness  and  high  frequency  of  the  source.  The 
waveforms  for  the  FE  solution  without  hourglass  corrections  show  marked  discrepan¬ 
cies  when  compared  to  either  the  corrected  FE  solutions  or  the  synthetics.  For  alter¬ 
nate  nodes  in  any  row  or  column  these  discrepancies  have  a  different  character.  In  one 
case  the  shape  of  the  waveform  shows  only  a  small  discrepancy  and  the  peak  ampli¬ 
tude  is  increased  by  more  than  10%.  In  the  other  case  the  entire  character  of  the 
waveform  changes  and  the  amplitude  decreases  by  at  least  20%.  Clearly  the  hourglass 
corrections  are  necessary  if  accurate  results  are  to  be  determined. 

The  second  example  used  a  100x100  node  FE  grid  whose  first  column  is  1500  km 
from  a  ten  kilometer  deep  line  source  in  a  layer  over  a  half  space.  The  layer  has  den¬ 
sity  2.7  g/cm3  and  SH  wave  velocity  3.5  km/s  and  the  half  space  has  density  3.2 
g/cm3  and  SH  wave  velocity  4.5  km/s.  The  time  step  size  and  the  vertical  and  hor¬ 
izontal  spacing  of  nodes  in  the  FE  grid  were  identical  to  those  used  in  the  previous 
example.  Seismograms  were  calculated  for  a  group  of  nodes  abng  the  surface  and 
down  depth  sections  at  twenty  five,  forty  and  sixty  five  nodes  from  the  grid  edge  using 
a  two  step  procedure.  First,  synthetic  seismograms  were  generated  at  a  horizontal 


distance  of  1500  km  from  the  source,  at  the  surface  and  at  half  kilometer  depth  inter¬ 
vals  to  a  depth  of  thirty  kilometers.  These  seismograms  were  used  to  constrain  the 
displacement  time  histories  of  the  leftmost  column  of  nodes  in  the  FE  grid.  It  will  be 
shown  that  applying  such  constraints  completely  specifies  the  distant  source.  Analytic 
synthetic  seismograms  were  also  generated  at  the  same  locations.  The  two  step 
seismograms,  which  have  been  propagated  through  a  FE  grid  for  thirty  to  thirty  five 
kilometers  between  1500  km  and  the  receiver  should  be  identical  to  the  analytic  syn¬ 
thetics  at  the  same  receivers  if  the  FE  calculation  is  stable.  Results  from  this  test  are 
shown  in  Figure  4.  The  pairs  of  two  step  seismograms  illustrated  show  how  the  intro¬ 
duction  of  the  hourglass  corrections  effects  the  waveform  and  amplitude  of  the  results. 
The  upper  trace  in  each  pair  includes  the  hourglass  restoring  force  the  lower  trace  does 
not.  At  the  scale  of  the  figure  the  corresponding  analytic  synthetics  are  indistinguish¬ 
able  from  the  two  step  results  with  hourglass  restoring  terms.  The  uppermost  number 
to  the  right  of  each  trace  is  the  ratio  of  the  peak  to  peak  amplitude  of  the  illustrated 
trace  to  that  of  a  synthetic  calculated  for  the  same  location.  The  lower  number  is  the 
same  type  of  ratio  using  an  RMS  amplitude  measure  discussed  in  detail  later.  These 
three  pairs  of  seismograms  were  recorded  at  the  free  surface  at  distances  of  30,  32.5, 
and  35  km  from  the  left  hand  grid  edge.  These  distances  illustrate  the  behavior  of  the 
amplitudes  at  adjacent  nodes  within  the  grid.  Some  intermediate  nodes  have  been 
omitted,  but  the  alternation  of  higher  and  lower  amplitudes  for  the  uncorrected  syn¬ 
thetics  is  seen  as  long  as  an  even  number  of  nodes  are  skipped  between  illustrated 
nodes.  Examining  the  illustrated  results  shows  that  adding  the  hourglass  restoring 
force  term  makes  a  small  but  perceptible  change  to  the  waveform.  This  change  pro¬ 
duces  a  slight  improvement  to  the  already  excellent  correspondence  between  the  FE 
and  analytic  waveforms.  However,  the  important  improvement  made  when  hourglass 
restoring  force  terms  are  included  is  in  the  amplitude  correspondence  between  the  FE 
and  synthetic  results.  When  no  hourglass  corrections  are  made  the  amplitudes  of  FE 


Figure  4:  Waveforms  for  hybrid  seismograms  calculated  in  a  FE  grid  removed 
from  the  source  by  a  distance  of  1500  km.  Seismograms  are  calculated  for  a  line 
source  in  a  layer  over  a  half-space.  Each  pair  of  seismograms  is  recorded  at  a  given 
node.  The  upper  trace  in  each  pair  includes  the  hourglass  restoring  force  terms.  The 
lower  trace  in  each  pair  does  not.  The  number  to  the  right  of  each  trace  are  ampli¬ 
tude  ratios.  The  amplitude  of  the  illustrated  seismogram  is  compared  to  that  of  a 
synthetic  calculated  for  the  same  location.  The  upper  number  is  a  ratio  for  peak  to 
peak  amplitudes,  the  lower  for  RMS  amplitudes.  Seismograms  are  of  duration  102 
seconds. 


seismograms  at  successive  nodes  along  any  row  or  down  any  column  alternate  between 
being  larger  and  smaller  than  those  of  the  corresponding  synthetics.  Introduction  of 
the  hourglass  restoring  force  terms  removes  this  oscillation  in  the  amplitude.  Clearly 
the  presence  of  the  oscillatory  term  caused  by  the  hourglass  instabilities  superimposed 
on  the  correct  solution  is  not  desirable.  Thus,  it  is  necessary  to  apply  the  hourglass 
restoring  force  correction  to  obtain  the  correct  solution. 

With  these  additions  the  SWIS  FE  code  produces  accurate  solutions.  However, 
the  calculations  are  extremely  time  consuming.  An  effort  was  made  to  increase  the 
speed  of  the  calculations  for  the  Cartesian  SH  option  used  for  the  calculations  in  this 
study.  First,  the  options  for  general  curvilinear  and  cylindrical  coordinates  were 
removed.  Some  small  subroutines  were  incorporated  into  higher  level  routines  to 
reduce  the  number  of  subroutine  calls.  The  system  was  defined  to  have  two  dimen¬ 
sions  and  one  degree  of  freedom  thus  removing  a  number  of  decision  statements  and 
blocks  of  code  necessary  to  the  implementation  of  other  options.  These  changes  also 
allowed  the  removal  of  a  number  of  do  loops  and  their  replacement  with  single  state¬ 
ments.  Other  modifications  were  made  where  incorporation  of  analytic  simplifications 
of  the  expression  being  evaluated  by  the  code  substantially  reduced  the  number  of 
necessary  calculations.  The  overall  effect  of  these  modifications  was  an  increase  of 
between  a  factor  of  four  and  and  factor  of  five  in  the  speed  of  execution  of  the 
modified  and  accelerated  SWIS  code  over  the  original  SWIS  code  for  case  of  the  propa¬ 
gation  of  SH  waves  in  Cartesian  coordinates. 

The  FE  and  propagator  matrix  methods  used  in  the  studies  have  now  been 
explained,  and  the  components  of  each  which  must  be  understood  to  explain  the 
Representation  Theorem  coupling  methods  have  been  described.  The  modifications 
made  to  the  SWIS  FE  code  have  increased  the  speed  of  the  code  by  more  than  a  fac¬ 
tor  of  five.  They  have  also  allowed  the  application  of  time  dependent  displacement 


constraints  to  selected  nodes  in  the  grid,  and  the  recording  of  time  slices,  or  snapshots 
of  the  motion  of  each  node  in  the  grid  at  a  given  time.  These  additional  capabilities 
are  utilized  to  implement  the  propagator  to  FE  coupling,  to  allow  the  use  of  larger  FE 
grids,  and  to  better  understand  the  propagation  of  disturbances  through  complex 
structures.  Next  it  is  necessary  to  understand  the  Representation  Theorem  (RT)  and 
how  it  can  be  applied  to  coupling  these  two  methods  to  produce  a  hybrid  method  util¬ 
izing  the  strengths  of  each. 

The  Representation  Theorem  and  Green’s  Functions 

The  basis  of  the  method  used  to  couple  the  propagator  matrix  calculation  to  the 
FE  calculation,  or  to  couple  the  FE  calculation  into  the  Propagator  matrix  calculation 
is  the  application  of  the  Representation  Theorem  on  the  boundary  between  the  regions 
in  which  each  method  is  used.  The  Representation  Theorem  (RT)  relates  the  displace¬ 
ment  at  a  point  in  a  volume  V  to  the  body  forces  fj  in  the  volume  V  and  to  the  dis¬ 
placements  u  and  the  tractions  Ts  on  the  surface  S  of  volume  V.  There  are  many 
equivalent  ways  of  expressing  the  RT,  for  example  (  Aki  and  Richards;  eq  2.41) 

-t-00 

up(*.t)  —  /  dr  J£f  fj(| ,  r)  Gip(X,t-r  ;  ?,0)  dV(|)  (26) 

4*00 

+  /  drff  jGip(2,t-r  ;  ?,0)  Ts  [  u(|,r),  n j  ]  -  cjjU  -^-Gkp(5,t-r  ;  J,0)  Uj(|  ,  r)  Dj  J  dS(0 

where  ^  the  location  of  a  point  of  the  surface  S  of  volume  V,  3?  the  location  of  a 
receiver  in  volume  V,  up(5,t)  is  the  p  component  of  the  displacement  at  time  t  at  the 
receiver  location,  ?,  t  is  the  time  at  which  the  observation  is  made,  r  is  the  source 
time,  Tj(?  ,  r)  is  the  boundary  condition  specifying  stress  as  a  function  of  source  time  r 
for  all  points  on  surface  S,  nj  is  the  jth  component  of  the  outward  unit  normal  to  the 
surface  S,  ut(? ,  r)  is  the  boundary  condition  specifying  displacement  as  a  function  of 
source  time  r  for  all  points  ^  on  the  surface  S,  and  Gjp  is  the  Green’s  function  which 
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represents  the  displacement  in  the  ith  direction  at  3?  at  time  t  due  to  a  unit  impulse 
applied  in  the  plh  direction  at  position  f  at  time  r  .  The  desired  Green’s  functions  are 
obtained  by  solving  the  wave  equation  below  in  volume  V,  subject  to  the  initial  condi- 

Q 

tions  that  GjP  and  -g— Gip  are  zero  for  t <r. 

6*  5(3?  -  ?)  5(t  -  r)  =  P~G-ip(S  ,  t  ;  ?  ,  r)  -  JL  { cijkl-£-Gkp(X ,  t  ;  1 ,  r)  J  (27) 

The  form  of  the  RT  given  as  equation  (26)  is  not  optimal  to  demonstrate  how  the 
coupling  of  the  two  methods  is  accomplished.  To  transform  it  to  a  more  manageable 
form  several  conditions  are  applied.  First,  it  is  assumed  that  no  body  forces  are 
present  within  volume  V.  Thus,  the  first  term  on  the  right  hand  side  of  equation  (26) 
vanishes.  Next,  the  medium  is  assumed  to  be  isotropic.  This  reduces  the  number  of 
independent  C;jkl  terms  to  nine  from  eighty  one.  It  allows  them  to  be  written  as 

Cjjkl  *“  *  *ij  *kL  +  f*  (*it  5jl  “  *U  ^jk)  (28) 

Using  these  values  of  cijkl  the  tractions  Tj  can  be  expressed  as 

Tj  =  Ty  n;  =  cijpq  Tj|-Up(?  .  r)  at  —  J  X  5jj  uu  +  M  (uu  +  ujti)  J  nj  (29) 
Substituting  (28),  and  (29)  into  (26)  a  more  useful  form  of  the  RT  is  found. 
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Gpirij  +  ^Gpij  Uj  ■+■  (Gpy  +  Gpjj)  Uj  J  nj  dS(£)  (30) 

Where  all  derivatives  are  with  respect  to  For  the  case  of  SH  waves  only,  (30)  can  be 
further  simplified  by  applying  the  conditions  governing  SH  motion  and  by  assuming 
all  quantities  are  constant  with  respect  to  x2,  that  is 
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e-1  dt 


«/  ,2  r2/a2 

r/fi  1  ~r  IP 


-J-f 

4ir^  J 


U./0U*  +  r«MB 


exp-^ - — — ^-d£ 

({!  +  r’)1'* 


where  r2  *«  (xj-^j)2  +  (xj-£3)2.  This  form  is  valid  for  a  whole  space,  and  must  be 


modified  to  satisfy  the  half-space  problem.  The  whole  space  problem  is  a  solution  to 


equation  (27)  subject  to  the  radiation  condition  that  rM-*0and  — - — — O  as  r— -oo. 

OX 


For  the  half-space  problem  rM  must  satisfy  an  additional  BC  as  well  as  satisfying 


equation  (27)  and  the  radiation  condition.  The  BC  must  satisfy  is  that  there  be  no 


stress,  <r~,  on  the  free  surface.  That  is,  the  relation 


0,  must  be  satisfied  on 


the  free  surface,  at  z=0.  This  is  equivalent  to  specifying  that  volume  V  is  a  half-space 


rather  than  a  whole  space.  Physically,  the  introduction  of  the  free  surface  boundary 


at  z=0  implies  that  energy  that  would  otherwise  continue  to  radiate  towards  infinity 


will  instead  be  reflected  back  toward  the  depth  of  the  receiver.  This  implies  that  the 


addition  of  a  term  of  the  same  form  as  the  whole  space  solution,  so  that  it  will  satisfy 


the  original  equation,  but  with  the  propagation  distance  r  equal  to  the  propagation 


path  length  of  a  surface  reflection  rather  than  a  direct  arrival,  will  satisfy  the  zero 


stress  BC  at  the  surface.  The  half-space  Green’s  function  for  a  SH  line  source  given 


below  is  based  on  this  reasoning  and  can  be  shown  to  satisfy  the  BC,  the  radiation 


condition,  and  the  initial  equation. 
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Uj  *  U3  —  0  (31) 

rU  —  TXS  ~  r13  *  r3I  “  0 
GjX  =»  Gja  ■■  G3j  a—  Gu  =«*  0 

_± _ 2__0 

dx2  d£2 

Also,  it  is  usual  to  assume  that  the  motion  and  stress  are  everywhere  zero  for  times 
less  than  f»0.  Applying  this  condition  and  the  conditions  given  in  (31)  to  (30)  and 
integrating  from  -oo  to  oo  over  x2  gives  the  form  of  RT  appropriate  for  a  2-D  SH 


line  source. 


0° 

«2(*.0  —  f*r  //*{r»uu  4-  r2S.ku8}nkdC(^e3r) 


where  C  is  the  curve  defined  by  the  intersection  of  the  surface  S  with  the  Xi~x3  plane, 
and  following  the  notation  of  de  Hoop  (de  Hoop  1958) 


rij  —  /  Gij(xi>x3,t)dxa 


Omitting  the  integration  over  x2  gives  a  SH  RT  appropriate  for  a  3-D  analogue  of  the 
procedure  described  below.  Further  discussions  in  this  chapter  relate  to  the  2-D  prob- 


The  form  of  the  line  source  Green’s  function  given  by  de  Hoop  (de  Hoop  1958)  is 


B 


where  r*  =  \f{*i-Zi)2  +  To  compare  this  form  for  the  Green’s  functions  to 

the  forms  of  the  displacements  used  in  the  rest  of  this  discussion  it  is  first  necessary  to 
transform  the  equation  from  the  Laplace  transform  domain  to  the  Fourier  transform 
domain,  and  to  express  the  modified  Bessel  function,  Kq,  in  terms  of  Hankel  functions 
or  of  exponentials.  Examining  the  second  expression  in  the  first  line  of  equation  (34) 
indicates  that  it  is  the  Laplace  transform  of  a  quantity,  f(r,z),  where  f(r,z)  is  the  first 
term  in  the  integrand.  Applying  the  change  of  variables,  £2  =  (vt)2-r2,  changes  the 
lower  limit  of  integration  to  zero  and  the  integration  variable  to  (.  Then,  noticing  the 
symmetry  of  the  resulting  integrand  allows  the  lower  limit  of  the  integral  to  be 
changed  to  -oo  if  a  factor  of  one  half  is  introduced.  This  gives  the  form  of  the  equa¬ 
tion  in  the  second  line  of  equation  (34).  Finally,  introducing  an  additional  change  of 
variables,  s  =»  iu,  changes  the  Laplace  transform  to  a  Fourier  transform.  The  result¬ 
ing  expression  for  the  Green’s  function  is 

r“  ”  2^  [  K«(,k»r*)  +  -  f  (M> 

Next,  it  is  useful  to  transform  the  modified  Bessel  functions,  K®,  into  Hankel  functions 
to  make  direct  comparison  with  the  displacement  solution  for  the  line  source  in  a  lay¬ 
ered  half-space  possible.  The  modified  Bessel  function  can  be  expressed  in  terms  of 
Hankel  functions  (Abramowitz  and  Stegun,  eq  9.6.4) 

KM  -  “  e-1"  H  « (ze-^2)  -~<*rg(z)<x  (-/) 

The  argument  of  the  modified  Bessel  function,  ik^r,  is  of  the  form  ix  where  x  is  a  posi¬ 
tive  real  number.  Thus,  arg(z)— |*,  and  the  expression  (37),  is  applicable.  Substitut¬ 
ing  equation  (37),  with  »/— 0,  into  equation  (36)  gives  the  first  expression  in  the  next 
group  of  equalities.  It  is  in  the  form  that  it  can  be  directly  compared  to  a  line  source 
displacement  expression  for  a  homogeneous  half-space  that  is  used  in  the  propagator 
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technique  to  develop  the  expression  for  the  layered  half-space.  This  expression  is 


*W.»)  -  Y  [Ho'^k^'  +  Ho'^kjr-'] 


Alternately,  the  equation  (36)  can  be  written  in  terms  of  integrals  of  exponentials  and 
then  be  compared  directly  to  the  solution  for  a  line  source  in  a  homogeneous  half¬ 
space.  The  exponential  form  of  the  equation  is  the  result  of  substituting  the  appropri¬ 
ately  transformed  version  of  the  integral  expression  in  equation  (34)  for  each  modified 
Bessel  function  in  equation  (36).  This  representation  of  the  Green’s  function  is  given 
as  the  second  equality  in  the  next  group  of  equations. 
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Comparing  the  second  equality  in  equation  (39)  to  equation  (25)  for  the  line 
source  in  a  homogeneous  half  space,  and  comparing  the  first  equality  in  (38)  to  the 
alternate  expression  for  the  line  source  in  the  homogeneous  half-space,  equation  (38), 
shows  that 

^(X,Z)“  (40) 

Thus,  the  displacement  Green’s  function  can  be  calculated  in  a  manner  identical  to 
the  displacement  solution,  in  either  a  layered  half-space  or  a  homogeneous  half  space. 

However,  one  must  be  careful  to  include  the  multiplicative  factor  of  when  deter¬ 
mining  the  Green’s  functions  for  use  in  the  Representation  Theorem  expressions. 
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The  form  of  the  Representation  Theorem  (RT)  given  in  equation  (32)  has  been 
used  to  determine  approximate  analytic  solutions  to  the  problems  of  propagation 
across  a  surface  step  (Knopoff  and  Hudson  1964),  and  of  propagation  in  a  layered  elas¬ 
tic  wedge  (Hudson  and  Knopoff  1964).  This  is  also  the  form  used  in  this  study  to 
develop  procedures  for  the  coupling  of  FE  and  Propagator  Matrix  solutions. 

Figures  5  and  6  show  geometries  useful  for  developing  and  explaining  the 
methods  to  be  used  for  coupling  FE  and  propagator  matrix  methods.  The  purpose  of 
developing  these  methods  is  to  allow  the  propagation  of  wavefields  from  a  source  to  a 
receiver  through  a  series  of  plane  layered  media  separated  by  complex  regions  that 
cannot  be  modeled  using  a  plane  layered  structure,  or  by  changes  in  the  layered  struc¬ 
ture.  For  example,  the  complex  regions  that  will  be  considered  in  chapter  2  are  con¬ 
tinent  to  ocean  and  ocean  to  continent  boundaries.  A  simple  geometry,  consisting  of  a 
layer  over  a  half-space,  is  used  as  a  test  case  to  show  how  the  methods  work,  and  how 
they  can  be  applied  in  both  this  simple  case  and  in  more  complicated  cases.  Figure  5 
can  be  used  to  illustrate  the  geometry  used  to  couple  energy  traveling  through  a  lay¬ 
ered  medium,  into  a  FE  grid  which  may  contain  structures  of  arbitrary  complexity. 
Figure  6  can  be  used  to  illustrate  the  coupling  from  such  a  FE  grid  into  a  layered 
structure  where  propagator  matrix  methods  are  used  to  transmit  the  energy  further. 
Each  of  these  figures  will  be  discussed  is  detail  later. 

Propagator  to  Finite  Element  Coupling 

In  a  layered  medium  the  wavefield  can  be  mathematically  constructed  at  any 
point  receiver  within  the  layered  half-space  using  the  propagator  matrix  technique  and 
an  appropriate  form  of  a  source  representation.  The  resulting  seismogram  at  any 
point  in  the  layered  half-space  will  include  not  only  a  direct  arrival  and  a  surface 
reflection  but  also  the  superposition  of  many  multiple  reflections  which  produce  the 
surface  waves  in  the  wavetrain.  Although  the  seismograms  are  produced  as  modal 
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sums,  in  the  limit  of  using  all  possible  modes,  they  are  equivalent  to  summing  over  all 
possible  multiple  reflections  in  terms  of  rays.  Thus,  they  can  be  conceptualized  in 
terms  of  rays.  To  simplify  the  graphical  representation  used  in  Figures  5  and  6  the 
seismogram  at  a  source  point  in  the  layered  half-space  is  represented  by  a  single  direct 
line  from  the  source  to  the  receiver.  This  line  also  implicitly  represents  all  the 
reflected  and  multiply  reflected  rays  in  that  seismogram  which  would  hopelessly 
clutter  the  diagrams  if  they  were  included. 

Figure  5  shows  a  graphical  representation  of  a  layer  over  a  half-space.  The 
source  is  represented  by  a  large  asterisk,  and  each  receiver  by  a  large  point.  The  ori¬ 
gin  of  the  coordinate  axes  is  also  shown  as  a  large  point.  The  coordinate  axes  Xj  =*  x, 
and  x2  »  z,  are  labeled  at  the  origin  of  the  Cartesian  coordinate  system  used  in  this 
discussion.  The  arrows  show  the  directions  in  which  the  values  of  each  coordinate 
increase.  The  distance  x  increases  to  the  right,  the  depth  z  increases  downwards,  and 
the  coordinate  y  **  x2  increases  out  of  the  page.  All  motion  occurs  in  the  y  direction 
since  this  is  an  SH  wave  problem.  The  free  surface  of  the  half-space  is  located  at  z=0, 
and  the  source  is  located  at  x*0  at  depth  h.  The  second  long  horizontal  line  is 
located  at  z—D,  where  D  is  the  thickness  of  the  layer  above  the  half-space.  The 
points  denoting  the  receivers,  also  illustrate  the  locations  of  the  nodes  in  two  columns 
of  a  FE  grid.  The  intervening  columns  of  nodes  are  not  illustrated  to  avoid  cluttering 
the  diagram.  The  short  horizontal  line  that  connects  the  lowermost  nodes  is  at 
z  *  Dre,  and  denotes  the  bottom  of  the  FE  grid.  The  FE  grid  is  assumed  to  extend 
off  the  page  toward  the  right.  The  vertical  line  connecting  the  leftmost  column  of 
receivers  shows  the  lefthand  edge  of  the  FE  grid.  For  the  sake  of  simplicity  and  clar¬ 
ity  the  structure  illustrated  within  the  FE  grid  is  the  same  layer  over  the  same  half¬ 
space  as  that  used  in  the  region  traversed  using  the  propagator  matrix  calculation. 


The  coupling  of  a  w&vefield  defined  by  the  discrete  sampling  of  seismograms 
along  the  leftmost  edge  of  a  FE  grid,  into  that  FE  grid  is  straightforward.  The 
seismograms,  for  a  specified  source  function,  are  generated  at  a  group  of  receivers 
equally  spaced  in  z  located  a  horizontal  distance  X  from  the  source.  The  depth  spac¬ 
ing  Az  between  the  receivers  is  the  node  spacing  in  the  FE  grid  into  which  the 
wavefield  is  to  be  coupled.  The  distance  X  is  the  minimum  separation  between  the 
source  and  the  lefthand  edge  of  the  grid.  The  solid  lines  connecting  the  source  to  the 
receivers  along  the  left  edge  of  the  grid  illustrates  these  seismograms.  They  are  calcu¬ 
lated  by  transmitting  the  appropriate  source  functions  through  the  layer  over  a  half¬ 
space  structure  from  the  source  to  the  leftmost  edge  of  the  FE  grid  by  using  the  prop¬ 
agator  matrix  method.  The  seismograms  are  then  applied  as  displacement  time  his¬ 
tory  constraints  on  the  leftmost  column  of  nodes  in  the  FE  gTid.  The  application  of 
this  type  of  constraints  to  a  column  of  nodes  in  the  FE  grid  completely  specifies  the 
subsequent  motion  at  all  points  in  that  FE  grid.  Thus,  a  hybrid  seismogram  can  be 
recorded  at  any  node  in  the  FE  grid  provided  a  FE  calculation  of  suitable  duration  is 
completed.  Consider,  as  an  example,  that  the  receiver  at  which  the  hybrid  and  ana¬ 
lytic  seismograms  are  to  be  compared  is  the  surface  receiver  at  a  distance  X2  from  the 
source.  The  minimum  distance  from  the  source  to  the  second  column  of  nodes  illus¬ 
trated  in  Figure  S  is  X2.  The  heavy  solid  line  in  Figure  5  represents  the  analytic 
seismogram,  calculated  using  only  the  propagator  matrix  method  and  the  appropriate 
source  functions.  The  hybrid  seismogram  is  represented  by  the  integration  over  all  the 
dashed  paths  connecting  the  edge  nodes  and  the  receiver.  Each  dashed  path 
represents  a  seismogram  for  which  the  appropriate  forcing  function  is  a  time  depen¬ 
dent  source. 

Figure  5  shows  the  layer  over  a  half-space  structure  extending  into  the  FE  grid. 
A  series  of  calculations  completed  using  such  a  homogeneous  structure  provides  a 


useful  test  of  the  procedures  used.  It  allows  the  comparison  of  the  hybrid  solutions, 
those  solutions  propagated  from  the  source  to  the  receiver  using  a  combination  of 
methods,  to  analytic  synthetics,  those  calculated  entirely  with  the  source  functions 
and  the  propagator  matrix  technique,  for  receivers  in  locations  identical  with  respect 
to  the  source.  Hybrid  and  analytic  synthetics  calculated  by  using  the  layer  over  a 
half-space  structure  in  both  portions  of  the  path  will  be  shown  later.  These  seismo¬ 
grams  will  be  used  to  demonstrate  the  validity  of  the  coupling  method  by  comparing 
the  hybrid  results  to  the  analytic  results.  Although,  in  this  example,  the  FE  grid  con¬ 
tains  the  same  layer  over  a  half-space  structure  as  that  used  for  the  propagator  matrix 
calculations,  any  other  arbitrary  structure  can  be  put  into  the  FE  grid  without  alter¬ 
ing  the  coupling  method.  However,  the  calculation  of  analytic  synthetics  to  which  the 
resulting  hybrid  synthetics  can  be  compared  may  be  difficult  or  impossible.  Thus, 
although  we  wish  to  address  the  effects  of  complex  structures  within  the  FE  portion  of 
the  path  the  simple  esse  where  the  FE  grid  contains  the  same  layered  structure  as 
used  for  the  propagator  matrix  calculation  will  be  discussed  here  as  it  is  the  best  test 
of  the  method. 

In  the  examples  discussed  in  later  sections  the  seismograms  used  as  forcing  func¬ 
tions  are  generated  using  more  than  one  type  of  source.  The  first  types  of  sources 
used  are  line  sources,  for  a  single  point  force  in  a  homogeneous  half-space  or  a  layered 
half-space.  For  these  cases  the  applied  forcing  functions  are  u7(x,y,z).  However,  the 
balance  of  the  sources  used  are  double  couple  point  sources.  The  propagator  matrix 
solutions  for  these  sources  are  of  the  form  7(r,4,z).  The  displacements  are  expressed  in 
cylindrical  coordinates  rather  than  in  the  Cartesian  coordinates  appropriate  to  the  RT 
integral  or  the  Cartesian  FE  representation.  Thus  it  must  be  demonstrated  that  these 
cylindrical  displacements  can  be  used  in  place  of  the  Cartesian  ones  expected  without 
adversely  affecting  the  results.  This  will  be  demonstrated  in  the  following  paragraphs. 


The  expressions  for  V(r,^,z)  for  the  dip  slip  and  strike  slip  faults,  equations  (6), 
are  in  cylindrical  coordinates  (r,^,z),  where  r  =  \/(x2+y2).  Displacement  time  histories 
of  a  series  of  locations  corresponding  to  the  first  column  of  nodes  in  a  FE  grid  a  dis¬ 
tance  rt  from  the  source  were  evaluated  using  these  expressions.  When  these  displace¬ 
ment  seismograms  are  used  as  input  forcing  functions  then  they  are  used  as  if  V(r,^,z) 
was  an  expression  for  ff,(x,y,z)  rather  than  for  TT^r.^.z).  No  transformation  from 
cylindrical  to  Cartesian  coordinates  is  performed.  The  displacement  field  defined  by 
the  imposition  of  displacement  time  history  constraints  down  the  FE  grid  edge  is  then 
propagated  an  additional  distance,  Ar,  through  the  FE  grid  to  give  the  hybrid  seismo¬ 
grams  that  approximate  IT#  at  the  receiver.  The  validity  of  using  V(r,^,z)  rather  than 
vTj(x,y,z)  will  now  be  demonstrated.  It  will  be  shown  that  for  ri»Ar  substituting 
u^  for  u;  gives  a  good  approximation  of  u*  at  the  receiver  despite  the  fact  that  the  FE 
method  produces  2-D  rather  than  3-D  propagation  effects. 

Consider  a  FE  grid  with  its  leftmost  edge  a  distance  from  a  source,  and  a 
receiver,  where  hybrid  and  analytic  synthetic  results  are  recorded,  a  distance  r2  from 
the  source.  Define  the  distance  propagated  within  the  FE  grid  as  Ar  =  r2  -  rr  Con¬ 
sider,  also,  the  expressions  (6)  for  the  displacements  from  a  point  dip  slip  or  strike  slip 
source.  Only  variations  in  the  r  coordinate  need  to  be  discussed  to  establish  the  vali¬ 
dity  of  this  procedure  so,  for  clarity,  (6a)  or  (6b)  will  be  reduced  to, 


<?Hi2>(kLr) 

dr 
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where  v  *  2  for  a  strike  slip  source  and  v  =  1  for  a  dip  slip  source.  The  derivative 
of  the  Hankel  function  can  be  expanded  in  terms  of  undifferentiated  Hankel  functions 
using  (Abramowitz  and  Stegun,  eq  9.1.29) 

•i-  H«HM  -  -  H«>,(M  +  -p-H  «'(kur) 
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(42) 


Each  resulting  Hankel  function  can  be  expanded  in  terms  of  the  asymptotic  expansion 
of  H„  for  large  r  (Abramowitz,  Stegun,  eq  9.2.4) 

Hj2)(kLr)  **  yj -£—■  e  (  V  -2ff<arg(kLr)<  jr  (43) 

Replacing  the  derivative  of  the  Hankel  function  in  equation  (41)  with  the  product  of 

kL  and  the  right  side  of  expression  (42),  gives  an  expression  for  jv(r>,z)  j  including 

only  undifferentiated  Hankel  functions.  Substituting  the  asymptotic  expansion,  equa¬ 
tion  (43),  for  each  of  these  Hankel  functions  gives  an  alternate  expression  for 

|v(r,^,z)J.  It  is 

{v(r>,z)J  —  kL“7  e  4  e_,k,/  (44) 

Since  we  are  considering  the  case  of  r  large,  kL»-^-,  and  the  second  term  in  the 

square  brackets  can  be  ignored,  jv(r2l^,r)J  can  be  expressed  in  terms  of  <t>  ,  z,  r,, 
and  Ar  as 


f  2kL  -4k,,.  -ik^r 
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(45) 
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Now,  for  a  line  source  in  a  half  space,  equation  (25),  or  for  a  line  source  in  a  vertically 
inhomogeneous  half-space,  the  modal  relation  analogous  to  equation  (45)  is, 


■j  !'J  "Aim "/  "JWV  v  v  vwv  v 
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u^(x2,z)  =  u#(xl,z)  *»  u#(xl,z)  e"*^  (46) 

Comparing  (46)  and  (45)  shows  that  both  expressions  have  the  same  form.  In  each 
case  the  displacement  at  rj  can  be  expressed  as  the  displacement  at  r]  multiplied  by  a 
propagation  factor.  A  FE  or  RT  calculation  will  give  the  same  propagation  factor  for 
each  mode  as  the  analytic  expression,  (46),  above.  If  the  displacements  at  rt, 

Vfo.^.z)  or  u,(x,z),  and  the  propagation  factors,  e“,kl'Ar  or  are  correlated  by 

considering  xl—rj  and  Ax=Ar,  then  V(ri,^,z),  the  3-D  solution,  will  be  given  by 
ur(x,z),  the  result  of  the  2-D  FE  or  RT  calculation,  so  long  as  the  source  is  many 
wavelengths,  k^r,  from  the  boundary  and  r  is  normal  to  the  FE  grid  edge  or  the  RT 

integration  surface.  1 
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Finite  Element  to  Propagator  Coupling 

The  coupling  of  FE  seismograms  into  a  layered  media  through  which  they  are 
transmitted  by  the  propagator  matrix  technique  and  the  application  of  the  RT,  is 
accomplished  by  direct  evaluation  of  the  RT  integral  as  given  in  equation  (32)  using 
propagator  matrix  generated  line  source  Green’s  functions.  The  geometry  of  the  prob¬ 
lem  is  illustrated  in  Figure  8.  The  geometry  will  be  explained,  then  the  coupling  pro¬ 
cedure  will  be  discussed. 

Figure  6  shows  the  geometry  used  to  evaluate  the  RT  integral,  equation  (32). 
Again,  the  layer  above  a  half-space  model  is  illustrated  in  both  the  FE  and  propagator 
matrix  regions  for  clarity.  In  practice  any  structure  can  be  used  in  the  FE  grid.  The 
two  long  horizontal  lines  show  the  free  surface  and  the  boundary  between  the  layer 
and  the  half-space.  The  solid  vertical  line  which  does  not  have  dots  superimposed 
upon  it  is  the  righthand  edge  of  the  FE  grid.  The  FE  grid  is  assumed  to  continue  off 
the  figure  to  the  left.  The  source  can  be  assumed  to  be  either  in  the  FE  grid  itself  or 
on  the  far  side  of  the  FE  grid.  The  short  horizontal  line  shows  the  bottom  of  the  FE 
grid.  The  x’s  shown  within  the  grid  and  on  its  bottom  boundary  represent  nodes  in 
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Figure  6:  Geometry  used  to  explain  the  coupling  of  FE  stress  and  displacement 
time  histories  into  propagator  matrix  calculations.  The  RT  integration  path,  C,  is 
shown  as  the  closed  curve  FGHF.  Finite  element  nodes  are  shown  as  open  circles.  The 
short  vertical  and  horizontal  lines  are  the  boundaries  of  the  FE  grid.  The  two  long 
horizontal  lines  show  the  free  surface  and  the  boundary  between  the  layer  and  the 
half-space.  The  filled  circles  superimposed  on  the  FG  section  of  the  integration  path 
are  element  centers  within  the  FE  grid  where  stress  and  displacement  time  histones 
are  recorded.  The  dot  within  the  integration  path  denotes  the  receiver.  The  lines  con¬ 
necting  the  element  centers  and  the  receiver  represent  the  line  source  Green's  func¬ 
tions. 


the  FE  grid.  The  two  illustrated  columns  of  nodes  define  the  corners  of  a  column  of 
elements.  The  center  of  each  of  these  elements  is  shown  as  a  dot.  The  integration 
path,  C,  for  the  RT  integral  is  illustrated  as  the  closed  curve  FGHF.  The  arrows  at  G 
and  H  indicate  that  the  integral  is  evaluated  on  the  contour,  C,  where  the  x  coordi¬ 
nate  at  G,  and  the  z  coordinate  at  H  both  tend  to  infinity.  An  example  of  a  possible 
receiver  inside  the  contour  C  is  illustrated  as  a  large  dot  outside  the  FE  grid  The 
lines  connecting  the  dots  showing  the  column  of  element  centers  within  the  FE  grid  to 
the  receiver  represent  the  line  source  Green’s  functions  for  propagation  from  a  source 
at  the  location  of  the  element  center  to  the  receiver. 

The  evaluation  of  the  RT  integral  on  the  contour,  C,  will  be  discussed  in  three 
parts.  First  the  quantities  used  in  the  evaluation  of  the  RT  integral  on  the  segment 
FG  of  the  contour,  C,  will  be  discussed.  Then,  the  procedure  used  to  estimate  the 
value  of  this  integral  will  be  explained.  Finally,  it  will  be  shown  that  the  contribu¬ 
tions  to  the  RT  integral  from  integration  along  segments  GH  and  HF  of  the  contour, 
C,  do  not  contribute  to  the  solution. 

Before  the  particulars  of  the  integration  procedure  can  be  discussed,  each  of  the 
quantities  in  equation  (32)  must  be  defined.  The  integration  surface  for  the  segment 
GH  is  a  line  with  the  x  coordinate  held  constant,  only  the  terms  generated  by  setting 
k»l  in  equation  (32)  need  to  be  considered.  Thus,  the  quantities  of  interest  are, 
H,  uy,  tJyrj,  and  ra>2.  The  displacement  and  stress  time  histories  at  the  element 
centers,  illustrated  as  dots  in  Figure  6,  are  recorded  during  the  FE  calculation.  The 
element  center  displacement  u,,  and  the  element  center  stress,  are  used  as  the 
u2  and  pu2>1  terms  respectively,  in  the  RT  integral  (32).  Line  source  Green’s  functions 
are  calculated  for  the  transmission  of  a  unit  line  displacement,  applied  at  each  of  the 
nodes  illustrated  as  dots  in  Figure  6  at  time  t,  to  the  receiver  point.  These  Green’s 
functions  are  calculated  using  the  propagator  matrix  method.  These  displacement 
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!  Green’s  functions  are  the  T terms  in  equation  (32).  The  evaluation  of  the  r2;2 

terms  is  analogous  to  the  calculation  of  the  displacement  Green’s  functions.  The 
details  of  this  calculation  will  be  discussed  later  when  the  evaluation  of  stress  seismo¬ 
grams  is  discussed.  Now,  all  the  necessary  time  dependent  quantities  in  the  RT 
integral  have  been  defined.  It  remains  to  note  that  the  rigidity,  n,  in  equation  (32)  is 
the  rigidity  at  the  depth  of  the  source. 

The  next  step  is  to  explain  how  the  RT  integral  is  numerically  integrated.  Con¬ 
sider  the  displacement  and  stress  seismograms,  and  the  displacement  and  stress 
Green’s  functions  as  time  series.  Each  time  series  gives  the  amplitude  of  a  displace¬ 
ment  or  a  stress  as  a  function  of  time.  The  time  spacing  between  successive  points  in 
each  time  series  is  the  time  step  duration,  At,  used  in  the  FE  calculation.  It  should  be 
remembered  that  the  terms  and  ra>kU2  >  >n  the  RT  integral,  equation  (32),  are 

convolutions.  Thus,  the  evaluation  of  the  integrand  at  each  integration  point  is  a  four 
step  process.  First,  each  of  the  displacement  and  stress  seismograms  is  Fourier 
transformed  into  the  frequency  domain.  Then,  the  products  of  the  Green’s  functions 
and  FE  results,  r^TT^i  aQd  are  calculated.  In  these  expressions  a  bar  denotes 

a  Fourier  transform.  Next,  the  resulting  products  are  inverse  Fourier  transformed 
back  into  the  time  domain.  Finally,  the  sum  of  the  two  convolutions  is  determined 
and  multiplied  by  the  rigidity.  Now,  the  value  of  the  integrand  at  each  element  center 
on  the  integration  surface  is  known.  Integration  along  the  segment  FG  or  the  contour, 
C,  is  approximated  by  taking  the  sum  over  the  values  of  the  integrand  at  each  ele¬ 
ment  center  along  the  subsegment  FFl  of  the  segment  FG  times  the  segment  length, 
and  applying  the  assumption  that  the  contributions  to  the  integral  from  the  remainder 
of  segment  FG  are  negligible.  The  assumption  that  no  significant  contribution  is 
made  by  integration  along  segment  FtG  can  be  justified  for  the  case  where  Ft  is 
chosen  so  that  the  earliest  possible  arrival  of  energy  from  a  source  at  depth  Fj  is  later 
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than  the  last  arrival  in  the  seismogram  being  considered.  For  the  remainder  of  this 
discussion  we  will  assume  that  the  depth  Ft  illustrated  in  Figure  6  satisfies  this  cri¬ 
terion.  Thus,  the  integration  of  the  RT  integral  on  the  segment  FG  of  the  contour,  C, 
is  expressed  as  a  sum  over  the  value  of  the  integrand  at  each  of  the  element  centers 
along  subsegment  FFt  times  the  segment  length.  This  is  equivalent  to  applying  a 
trapezoidal  rule  numerical  integration  to  the  equation  (32)  along  segment  FFt. 

Next,  the  integration  over  the  remaining  two  segments  of  the  contour,  C,  will  be 
shown  to  give  no  contribution  to  the  RT  integral.  The  segment  of  the  contour  FH, 
along  the  free  surface  will  be  considered  first.  The  integration  surface  for  the  segment 
FH  of  the  contour,  C,  is  a  line  with  the  z  coordinate  held  constant.  Thus,  only  terms 
generated  by  setting  k=3  in  equation  (32)  need  to  be  considered.  The  free  surface 
boundary  condition  applied  on  this  surface  states  that  the  is  zero  stress  on  this  sur¬ 
face.  Thus,  at  z=0  the  stress,  /iujj,  is  zero.  Since  the  Green’s  function,  r^,  used  in 
all  the  calculations  also  satisfies  the  free  surface  BC,  the  Green’s  function  stress, 
®  also  zero  at  z— *0.  Therefore,  one  term  in  each  product  in  the  integrand  of 
equation  (32)  is  zero,  causing  the  value  of  the  integral  along  this  portion  of  the  con¬ 
tour,  C,  to  be  zero.  Next,  the  segment  GH  of  the  contour,  C,  will  be  considered.  As 
r— *oo  the  displacements  r«  and  u2  must  approach  zero.  The  radiation  BC  used  in 
deriving  the  expressions  for  these  values  imposes  such  a  condition.  Substituting  the 
asymptotic  expansion  for  H„  into  equation  (38)  and  (39)  and  taking  the  limit  as  r-*oo 
also  verifies  this  statement.  Again,  one  component  of  each  product  in  the  integrand  of 
equation  (32)  is  zero.  Clearly,  this  makes  the  integrand  zero  and  verifies  that  the  sec¬ 
tion  GH  of  contour,  C,  makes  no  contribution  to  the  RT  integral. 

In  summary,  the  mechanics  of  the  methods  used  to  couple  propagator  matrix 
solutions  into  FC  calculations  and  to  couple  FE  results  into  the  RT  integral  using  line 
source  Green’s  functions  generated  with  the  propagator  matrix  method  are  described 
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briefly .  To  couple  propagator  matrix  solutions  into  a  FE  calculation  three  basic  steps 
are  necessary. 

(1)  Generate  seismograms  at  an  appropriate  set  of  receivers  using  the  desired  source 
functions  and  the  propagator  matrix  technique.  An  appropriate  set  of  receivers 
consists  of  receivers  at  a  fixed  horizontal  distance  from  the  source,  at  depths 
corresponding  to  the  nodes  in  the  first  column  of  the  FE  grid  to  be  used  in  the 
FE  calculation. 

(2)  Use  this  set  of  seismograms  as  displacement  time  history  constraints  on  the  left¬ 
most  column  of  nodes  in  that  FE  grid. 

(3)  Record  hybrid  seismograms  at  the  desired  locations  within  the  FE  grid  as  the 
constrained  FE  calculation  proceeds. 

To  couple  results  from  a  FE  calculation  back  into  a  simple  layered  structure  across 
which  they  can  be  transmitted  using  propagator  matrix  Green’s  functions  and  the  RT 
integral  five  steps  are  needed. 

(1)  Record  element  center  displacement  and  stress  time  histories  from  the  FE  calcu¬ 
lation. 

(2)  Calculate  line  source  displacement  and  stress  Green’s  functions  that  will  transmit 
a  disturbance  from  each  element  center  to  the  desired  receiver. 

(3)  Fourier  transform  the  displacement  and  stress  time  histories  and  the  displace¬ 
ment  and  stress  line  source  Green’s  functions  into  the  frequency  domain.  Per¬ 
form  the  multiplication  of  the  displacement  time  history  with  the  stress  Green’s 
function  and  the  stress  time  histories  with  the  displacement  Green’s  functions. 
Inverse  Fourier  transform  the  products  back  into  the  time  domain. 

(4)  Add  the  appropriate  pairs  of  products  and  multiply  by  the  rigidity  to  give  the 
values  of  the  integrand  at  each  of  the  element  centers. 


(5)  Sum  the  values  of  the  integrands  multiplied  by  their  respective  interval  lengths 

along  the  integration  surface  to  give  the  hybrid  seismogram. 

Tests  of  Accuracy:  SH  Pulse  Input 

The  firsthand  the  most  basic,  test  of  the  validity  of  the  RT  formulation  is  the 
coupling  of  two  successive  FE  calculations.  This  type  coupling  can  in  itself  be  used  to 
advantage  under  some  conditions,  and  will  be  referred  to  as  grid  extension  in  later  dis¬ 
cussions.  For  example,  propagation  through  successive  FE  grids  can  be  useful  in  some 
cases  to  deal  with  unwanted  boundary  reflections.  As  will  be  demonstrated  later  the 
transmitting  BC  used  is  angularly  dependent,  and  removes  little  of  the  reflection  at 
grazing  or  near  grazing  incidence.  For  FE  calculations  with  large  aspect  ratios,  that  is 
nx  >  >  nz,  wide  angle  reflections  from  the  grid  bottom  can  contaminate  the  results. 
When  seismograms  of  short  duration  are  required  these  wide  angle  reflections  can  be 
avoided  by  using  a  series  of  small  aspect  ratio  grids.  For  example,  if  a  series  of  two 
small  aspect  ratio  FE  calculations  are  used  to  model  a  long  aspect  ratio  problem  then 
it  should  be  possible  to  shift  at  least  half  the  width  of  the  grid  without  worrying 
about  wide  angle  reflections.  In  that  case  two  small  grid  runs  of  the  FE  would  be  fas¬ 
ter  than  one  run  with  an  expanded  grid.  Thus,  the  grid  extension  would  be 
worthwhile. 

In  Figure  7  a  pair  of  fifty  by  fifty  node  grids  are  illustrated.  A  line  force  source, 
which  is  a  point  force  when  it  is  projected  onto  the  x-z  plane  defined  by  the  FE  grid, 
is  applied  to  the  first  node  in  the  nlk  row  of  nodes  in  the  first,  or  upper  grid.  It  is 
shown  in  the  figure  as  an  crossed  circle  on  the  left  edge  of  the  upper  grid.  The  wide 
vertical  line  within  this  grid  represent  the  sixteenth  column  of  elements,  where  inter¬ 
mediate  seismograms  are  recorded  for  later  use  as  forcing  functions  input  to  the 
second  step  of  the  procedure.  The  narrow  vertical  line  in  the  upper  grid  represents 
the  thirty  first  column  of  nodes  where  direct  seismograms  representing  results  for  a 
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Figure  7:  Grids  used  to  test  coupling  of  two  successive  FE  calculations.  Both 
grids  are  illustrated  to  scale.  The  horizontal  lines  ending  in  arrows  are  labeled  with 
the  number  of  element  widths  they  span.  The  vertical  lines  within  the  grids  show  the 
locations  of  the  depth  sections  where  results  were  recorded.  The  crossed  circle  on  the 
left  side  of  the  first  grid  represents  the  line  source.  Both  grids  have  dimensions  50x50. 
The  horizontal  shift  of  the  second  grid  with  respect  to  the  first  grid  is  15  element 


single  FE  calculation  in  a  fifty  by  fifty  grid  are  recorded.  Figure  7  also  shows  the  spa¬ 
tial  relationship  of  the  grid  used  for  the  second  step  of  the  procedure  to  the  grid  used 
to  calculate  the  direct  results  and  the  forcing  functions.  The  first  column  of  nodes  in 
the  second  grid  is  shifted  to  the  right  to  coincide  with  the  location  of  the  sixteenth 
column  of  nodes  in  the  first  grid.  The  vertical  position  of  the  free  surface  in  each  grid 
is  identical.  For  the  second  step  the  intermediate  seismograms  recorded  in  step  one 
are  used  as  forcing  functions  on  the  first  column  of  nodes  in  the  second  grid.  The  sin¬ 
gle  vertical  line  within  the  second  grid  represents  column  sixteen,  where  ’two  step’ 
seismograms  are  recorded.  The  total  propagation  distance  for  the  'two  step’  pro¬ 
cedure,  summed  over  both  steps,  is  equal  to  the  distance  propagated  in  the  direct  cal¬ 
culation.  The  application  of  the  displacement  time  histories  used  ss  forcing  functions 
to  the  left  hand  edge,  or  the  first  column  of  nodes,  of  the  second  grid  demonstrates 
that  specifying  the  the  displacements  as  a  function  of  time  due  to  a  given  source,  for 
each  of  the  nodes  in  the  leftmost  column  of  a  FE  grid  constrains  the  motions  of  all 
points  in  the  grid  to  be  consistent  with  that  source.  When  any  column  other  than 
that  closest  to  the  source,  in  this  case  the  leftmost  column,  is  used  as  the  input 
column,  all  nodes  more  distant  from  the  source  than  the  input  column  are  consistent 
with  the  source.  The  nodes  source  ward  of  the  input  column  form  a  mirror  image  of 
the  distant  nodes  with  the  symmetry  line  being  the  input  column.  Any  energy  incident 
upon  the  boundary  formed  by  the  input  column  is  reflected  as  it  would  be  from  a  rigid 
boundary. 

Figure  8  shows  a  comparison  of  the  seismograms  generated  by  a  direct  calculi 
tion  and  by  the  ’two  step’  procedure.  The  top  seismogram  in  each  column  is  recorded 
at  the  surface.  The  vertical  spacing  between  the  nodes  at  which  the  remaining  seismo¬ 
grams  are  recorded  is  0.4  km  or  two  element  widths,  depth  increases  as  one  moves 
down  each  column.  The  first  column  shows  the  direct  results  recorded  at  column 


Figure  8:  Sample  waveforms  recorded  along  the  depth  section  thirty  element 
widths  from  the  source  that  is  illustrated  in  figure  1.2d.  The  first  column  shows  the 
direct  seismograms,  those  recorded  in  the  grid  containing  the  source.  The  second 
column  shows  the  two  step  results  recorded  at  the  same  distance  from  the  source  but 
in  the  grid  shifted  away  from  the  source.  The  vertical  spacing  between  nodes  where 
seismograms  are  illustrated  is  2  element  widths.  The  uppermost  pair  of  seismograms 
are  recorded  at  the  surface.  Depth  increases  as  one  moves  down  each  column. 


thirty  two  in  the  first  grid.  The  second  column  shows  the  corresponding  ’two  step’ 
results  recorded  at  column  sixteen  of  the  second  grid.  The  initial  waveforms  in  each 
column  are  almost  identical,  the  only  discernible  difference  being  a  second  pulse  result¬ 
ing  from  reflection  from  the  right  hand  end  of  the  grid  which  is  seen  for  the  direct  cal¬ 
culation  only.  The  reflection  is  not  seen  for  the  'two  step’  procedure  because  the  right 
hand  end  of  the  grid  is  fifteen  element  widths  farther  to  the  right,  moving  the 
reflection  to  a  time  later  than  the  end  of  the  seismogram.  The  amplitudes,  measured 
as  the  height  of  the  first  peak,  are  equivalent.  The  amplitude  variations  between 
methods  are  less  the  0.1%.  These  results  demonstrate  the  validity  of  coupling  two 
successive  FE  calculations  by  applying  forcing  functions  as  BC’s  on  the  edge  of  the 
second  FE  grid. 

Figure  9  shows  the  grid  geometry  used  to  illustrate  how  forcing  functions  can  be 
used  to  couple  energy  from  a  distant  source  into  a  FE  grid.  Seismograms  for  a  line 
source  2.5  km  below  the  free  surface  in  a  half-space,  at  a  horizontal  distance  of  fifteen 
kilometers  from  the  edge  of  the  grid,  were  determined  using  a  numerical  evaluation  of 
the  analytic  solution  for  the  direct  ray  and  the  ray  reflected  from  the  free  surface  (eq 
25).  The  seismograms  form  a  complete  set  of  forcing  functions  including  a  seismogram 
for  each  node  in  the  first  column  of  the  grid.  Analytic  seismograms  determined  by 
evaluation  of  the  expressions  for  a  line  source  in  a  half-space  (eq  25)  were  calculated  at 
horizontal  distances  from  the  source  corresponding  to  rows  sixteen  and  thirty-one  of 
the  grid.  Rows  sixteen,  thirty  one,  and  forty  six  are  shown  in  Figure  9  as  the  three 
vertical  lines  within  the  grid.  Using  the  forcing  functions  as  input  to  a  FE  calculation 
hybrid  solutions  were  calculated  for  all  nodes  in  columns  sixteen  and  thirty-one. 
Thus,  the  analytic  solutions  can  be  used  as  a  reference  against  which  the  hybrid  solu¬ 
tions  can  be  compared. 


Figure  9:  Geometry  vied  to  show  the  validity  of  coupling  energy  from  a  distant 
source  into  a  FE  calculation.  The  grid  illustrated  is  50x50  nodes,  and  is  shown 
scale.  The  source  is  illustrated  as  the  filled  circle  outside  the  grid.  The  source 
grid  distance  is  not  to  scale.  The  vertical  lines  in  the  grid  show  the  depth  sections 
where  seismograms  were  recorded.  These  seismograms  are  illustrated  in  Figure  10 
The  horizontal  lines  ending  in  arrows  are  labeled  with  their  lengths. 
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Figure  10  shows  the  results  of  the  test  of  the  validity  of  the  RT  coupling  of 
energy  from  a  distant  source  into  a  FE  calculation.  The  first  two  columns  show  the 
results  for  column  sixteen  of  the  grid,  which  is  three  kilometers  into  the  grid.  The 
second  two  columns  show  analogous  results  at  column  thirty  two,  six  kilometers  into 
the  grid.  The  first  column  in  each  group  shows  analytic  seismograms  for  reference. 
The  second  column  in  each  group  shows  the  hybrid  seismograms  at  the  corresponding 
locations.  The  first  seismogram  in  each  column  is  recorded  at  the  free  surface.  Suc¬ 
cessive  seismograms  moving  down  each  column  are  recorded  with  one  kilometer  depth 
spacing.  The  waveforms  of  the  analytic  and  the  hybrid  results  agree  extremely  well. 
However,  the  ’hybrid’  results  for  column  thirty  two  include  a  reflection  from  the  end 
boundary  of  the  grid,  not  present  in  the  analytic  calculations.  The  illustrated  seismo¬ 
grams  are  scaled  so  that  the  maximum  peak  has  unit  height.  Thus,  the  apparent 
amplitude  differences  in  the  initial  peaks  at  six  kilometers  is  an  artifact  of  the  scaling. 
The  correspondence  in  amplitude  between  the  two  methods  is  excellent.  At  three 
kilometers  the  agreement  deteriorates  below  the  fortieth  of  fifty  grid  points  from 
better  than  0.5%  to  as  much  as  10-20%.  At  six  kilometers  the  departure  from  the 
correct  amplitudes  begins  at  a  shallower  depth,  after  about  thirty  five  grid  points.  The 
variations  of  the  quality  of  fit  with  distance  from  the  grid  edge  and  with  depth  from 
the  free  surface,  can  easily  be  explained.  The  truncation  of  the  RT  integration  at  the 
bottom  of  the  grid  causes  artifacts  which  are  largest  as  the  depth  of  truncation  is 
approached  or  surpassed.  To  be  more  specific,  it  is  useful  to  extend  the  test  discussed 
above  to  more  clearly  delineate  the  nature  of  these  >  ^nations,  and  the  restrictions 
they  place  on  the  use  of  this  method  of  coupling  a  distant  source  into  a  FE  calculi 
tion.  A  more  extensive  test  will  be  discussed  below. 

To  minimize  the  calculation  needed  to  obtain  the  forcing  functions,  it  is  useful  to 
investigate  the  effect  of  reducing  the  number  of  forcing  functions  used,  that  is  limiting 
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Figure  10:  Demonstration  of  the  coupling  of  a  distance  source  into  a  FC  calcula¬ 
tion.  Two  pain  of  columns  showing  seismograms  recorded  at  the  depth  sections  three 
and  six  kilometen  from  the  grid  edge.  The  first  column  in  each  pair  shows  the  results 
of  the  FE  step  of  the  hybrid  calculation,  the  second  column  shows  the  half-space  syn¬ 
thetics  at  the  same  locations.  The  fint  row  of  seismograms  are  recorded  at  the  sur¬ 
face,  successive  rows  are  separated  by  one  kilometer  in  depth.  Peak  amplitudes  are 
normalized  to  one. 
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the  depth  extent  over  which  forcing  functions  are  applied  to  the  left  hand  grid  edge. 
The  exact  solution  of  the  RT  coupling  problem  requires  integration  to  infinite  depth. 
However,  energy  arriving  after  the  time  of  the  last  sample  in  the  modeled  seismogram 
can  be  ignored.  Thus,  the  maximum  depth  which  needs  to  be  considered  is  twice  the 
depth  of  the  deepest  receiver  plus  the  distance  over  which  the  travel  time  is  equal  to 
the  seismogram  duration.  It  is  clear  that  prohibitively  large  FC  grids  would  be 
required  to  produce  records  of  reasonable  length  or  at  large  distances.  However,  solu¬ 
tion  of  acceptable  accuracy  can  be  produced  by  integrating  over  a  comparatively  small 
number  of  forcing  functions.  The  previous  example,  of  coupling  energy  from  a  dis¬ 
tance  source  into  a  FC  grid  illustrates  that  this  is  the  case.  For  a  fifty  by  fifty  grid, 
discrepancies  of  less  than  0.3%  for  the  uppermost  forty  grid  points  are  observed.  A 
complete  solution  at  the  depth  of  the  fortieth  row  of  receivers  would  require  nz  to  be 
larger  than  one  hundred  fifty,  three  times  the  size  used  in  the  example.  This  observa¬ 
tion  implies  that  fewer  forcing  functions  are  needed  if  some  differences  at  depth,  near 
the  bottom  of  the  grid,  are  permissible.  An  experiment  to  investigate  the  effects  of 
applying  forcing  functions  to  only  the  top  n  nodes  in  the  input  column  has  been  per¬ 
formed.  A  geometry  and  procedure  analogous  to  that  used  to  demonstrate  the  vali¬ 
dity  of  the  RT  coupling  of  a  distant  source  into  a  FE  calculation  was  used. 

The  geometry  of  the  grids  used  to  investigate  the  effects  of  changing  the  number 
of  the  applied  forcing  functions  and  thus,  the  depth  where  the  deepest  forcing  function 
is  applied,  on  the  accuracy  of  the  resulting  seismograms  is  shown  in  Figure  11.  A  line 
source  embedded  in  a  half-space  at  a  depth  of  2.5  km  was  used  and  is  shown  as  a 
large  solid  dot  on  the  vertical  dotted  line  denoting  the  x  coordinate  origin.  The  hor¬ 
izontal  distance  from  the  source  to  the  left  hand  edge  of  the  grid  is  fifty  kilometers. 
The  first  column  of  nodes  in  the  grids  used  are  shown  as  a  pair  of  vertical  solid  lines. 
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Figure  11:  Grids  used  to  test  the  effect  of  varying  the  number  of  forcing  func¬ 
tions.  Double  vertical  lines  show  the  left  edge  of  the  grids.  Horizontal  lines  crossing 
these  vertical  lines  show  the  bottommost  forcing  function  for  each  test.  Tests  A15, 
A25,  A35,  and  A50  are  completed  in  the  50x50  grid  whose  extent  is  marked  by  the 
dotted  boundaries.  Tests  B30,  B50,  B70,  and  BlOO  are  done  in  the  large  grid  whose 
bottom  is  shown  as  a  horizontal  line  ending  in  a  solid  arrow.  The  solid  vertical  lines 
within  the  grid  show  the  depth  sections  where  results  are  recorded.  Solid  horizontal 
lines  ending  in  arrows  are  labeled  with  their  lengths. 


-116- 


synthetic  seismograms  are  compared  with  the  hybrid  seismograms  from  the  FE  calcu¬ 
lations.  They  are  at  distances  of  53.5  km  and  57.0  km  from  the  source,  corresponding 
to  columns  sixteen  and  thirty  one  in  the  FE  grids.  Two  grids  are  used  for  these  tests. 
The  dotted  lines  within  the  large  100x100  grid  are  the  bottom  and  right  hand  edges  of 
the  smaller  fifty  by  fifty  grid  used  for  the  first  four  tests.  The  solid  horizontal  lines 
ending  in  solid  arrows  show  the  free  surface  and  the  bottom  edge  of  the  large  FE  grid 
used  for  the  remaining  tests.  For  both  FE  grids  hybrid  seismograms  were  calculated 
for  four  cases.  Each  of  these  cases  used  a  different  number  of  forcing  functions.  For 
the  small  FE  grid  fifteen,  twenty  five,  thirty  five,  or  fifty  forcing  functions  were  used. 
For  the  large  FE  grid  thirty,  fifty,  seventy,  or  one  hundred  forcing  functions  were 
used.  The  vertical  extents  of  the  input  forcing  functions  are  illustrated  to  scale  in  Fig¬ 
ure  11.  The  locations  of  the  bottommost  input  forcing  function  for  the  small  grid 
tests  are  shown  by  the  horizontal  lines  labeled  A15,  A25,  A35,  and  A50  which  cross 
the  first  column  of  elements  in  the  illustrated  grid.  For  the  large  grid  the  bottommost 
forcing  functions  are  shown  as  the  lines  labeled  B30,  B50,  B70,  and  the  bottom  of  the 
large  grid,  B100.  Results  from  each  of  these  eight  tests  were  examined  to  determine 
the  effects  of  grid  size  and  vertical  extent  of  the  input  forcing  functions. 

Figure  12  shows  examples  of  the  seismograms  observed  at  column  sixteen  in  the 
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four  small  grid  tests.  The  surface  seismograms  are  illustrated  in  the  first  row  of  the 
figure.  Each  row  is  recorded  at  nodes  about  0.9  km  deeper,  that  is  four  element 
widths  deeper  than  the  previous  row.  The  column  of  numbers  in  the  middle  of  the 
figure  give  the  node  row  at  which  that  row  of  seismograms  was.  recorded.  The  first 
column  shows  results  from  test  A50.  These  waveforms  are  very  similar  to  the  analytic 
seismograms  recorded  for  the  same  locations.  The  agreement  is  best  at  the  surface,  at 
depths  below  the  thirty  fifth  node  the  decay  of  the  pulse  becomes  less  rapid  than  that 
seen  on  the  analytic  synthetic  and  the  maximum  amplitude  of  the  pulse  increases  with 


Figure  12:  Seismogram*  for  four  cases  using  different  numbers  of  forcing  func¬ 
tions.  Each  column  shows  seismograms  on  a  depth  section  at  53.5  km  from  the 
source.  The  first  row  is  recorded  at  the  surface.  Vertical  spacing  between  rows  is 
about  0.9  km.  The  first  column  uses  50  forcing  functions,  the  second  35,  the  third  25, 
and  the  fourth  15.  The  column  of  numbers  in  the  center  of  the  figure  gives  the  node 
row  number  for  each  row  of  seismograms. 


respect  to  the  analytic  synthetic.  For  this  test  the  difference  in  waveform  shape  is 
small.  The  second  column  shows  results  from  test  A35.  A  second  pulse  not  visible  for 
the  A50  test  is  present.  This  pulse  is  a  reflection  from  the  grid  bottom  of  energy  not 
canceled  by  the  disturbances  input  in  the  bottommost  fifteen  forcing  functions  in  test 
A50.  The  pulse  is  most  prominent  for  traces  recorded  at  depths  below  the  bottom¬ 
most  forcing  function  applied  in  this  test.  This  is  predictable  since  it  could  not  be 
expected  that  results  for  receivers  below  the  input  nodes  would  be  reliable.  If  the 
reflection  can  be  removed,  as  with  an  absorbing  BC  then  the  correspondence  of  A35 
results  to  direct  synthetics  is  good  down  to  node  row  thirty  or  thirty  five.  The  third 
column  shows  results  of  test  A25.  Again  a  reflected  pulse  contaminates  the  deeper 
receivers.  Results  compare  well  with  the  analytic  synthetics  down  to  about  node  row 
twenty.  The  fourth  column  shows  results  of  test  A15.  Here,  the  same  trend  already 
observed  continues.  The  correspondence  between  A15  results  and  analytic  synthetics 
is  good  down  to  about  node  row  ten.  It  appears  that  acceptable  waveforms  can  be 
calculated  for  depths  such  that  at  least  five  to  ten  input  forcing  functions  are  applied 
at  rows  below  the  receiver.  Similar  results  are  observed  for  these  tests  at  seven  kilom¬ 
eters  from  the  grid  edge,  and  for  the  large  grid  tests  at  both  distances.  Some  indica¬ 
tion  that  calculation  of  accurate  seismograms  at  longer  distances  may  require  a  larger 
vertical  extent  of  input  forcing  functions  has  been  observed. 

Figure  13  gives  a  more  complete  analysis  of  the  amplitudes  of  the  seismograms 
recorded  in  these  eight  tests  in  relation  to  the  amplitudes  of  analytic  synthetics.  At 
both  distances  the  maximum  amplitudes  of  the  seismograms  from  each  test  were  tabu¬ 
lated.  These  amplitudes  were  used  to  calculate  amplitude  ratios.  The  amplitude 
ratio,  A,  is  defined  as 
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Figure  13:  Relative  amplitudes  of  synthetic  and  two  step  seismograms  at  the 
same  locations,  as  a  function  of  the  number  of  input  forcing  functions.  The  left  hand 
graph  shows  results  recorded  at  a  depth  section  3.5  km  from  the  grid  edge,  the  right 
hand  graph  results  recorded  at  7  km.  Heavy  lines  show  tests  for  which  all  node  rows 
in  the  grid  have  forcing  functions  applied.  Lighter  solid  lines  show  results  of  tests  in 
the  50x50  grid.  Dotted  line  show  result  of  tests  in  the  100x100  grid.  The  number  in 
the  label  above  each  line  indicated  the  number  of  forcing  functions  used  in  that  test. 


These  amplitude  ratios  were  then  plotted  as  a  function  of  depth,  measured  as  a 
number  of  node  rows.  The  resulting  plots  are  shown  in  Figure  13.  The  left  hand  plot 
is  for  a  distance  of  3.5  Icm,  and  the  right  hand  plot  for  a  distance  of  seven  kilometers. 
The  heavy  lines  in  each  plot  indicate  the  amplitude  ratios  for  the  two  tests  in  which 
all  nodes  in  column  one  of  the  grid  have  forcing  functions  applied  to  them,  A50  and 
B100.  These  tests  are  used  as  a  reference  against  which  the  other  tests  can  be  com¬ 
pared.  The  dotted  lines  labeled  A15,  A25,  and  A35  show  where  the  amplitude  ratios 
for  the  three  tests  in  the  small  grid,  using  partial  sets  of  forcing  function  input,  depart 
from  those  of  test  A50.  The  amplitude  ratios  agree  well  with  those  of  the  A50  test  for 
the  near  surface  nodes  in  each  test.  The  depth  at  which  they  depart  from  the  A50 
results  increases  as  the  number  of  input  forcing  functions  increases.  As  expected  the 
amplitude  correspondence  is  better  when  more  complete  input  is  used.  Examining 
these  results  at  a  distance  of  3.5  km  shows  that  if  results  for  n  rows  are  desired,  input 
forcing  functions  should  be  applied  to  n+5  nodes.  At  seven  kilometers  similar  results 
are  observed  but  the  departure  from  acceptable  amplitude  ratios  occurs  four  or  five 
nodes  shallower.  Thus,  as  distance  increases  a  larger  depth  extent  of  input  forcing 
functions  is  needed  to  produce  amplitude  correspondence  within  acceptable  bounds. 
The  results  for  the  tests  using  the  larger  grid  are  show  in  Figure  13  as  the  light  lines 
labeled  B30  and  B50.  The  results  show  the  same  trends  as  the  results  for  the  small 
grid  indicating  that  the  artifacts  in  the  tests  using  truncated  input  sets  are  due  to 
interaction  with  the  grid  bottom  of  incompletely  integrated  signals  in  the  seismo¬ 
grams. 

If  one  wishes  the  analytic  solutions  to  correspond  well  to  the  FE  solutions  the  fre¬ 
quency  content  of  the  forcing  functions  used  to  drive  the  FE  calculation  becomes 
important.  To  investigate  the  effects  of  the  frequency  content  of  the  input,  sets  of 
forcing  functions  for  a  source  distant  from  the  FE  grid- were  calculated  for  two  cases. 


Each  of  these  sets  consists  of  seismograms  calculated  along  a  depth  section  fifty  kilom¬ 
eters  from  the  source,  with  a  vertical  node  spacing  of  0.2  km.  Both  sets  of  seismo¬ 
grams  were  determined  for  a  line  source,  with  a  triangular  time  function,  in  a  half¬ 
space  (eq  25).  The  first  case  used  a  time  function  with  a  duration  of  .15  seconds,  the 
second  used  a  time  function  with  a  duration  of  .45  seconds. 

The  results  of  these  tests  are  illustrated  in  Figures  14  15,  and  16.  The  first  two 
of  these  figures  show  the  seismograms  recorded  on  a  depth  section  53.5  km  from  the 
source.  The  vertical  spacing  between  successive  rows  of  seismograms  is  about  one 
kilometer.  The  amplitudes  in  all  traces  in  these  figures  have  been  normalized  such 
that  the  peak  amplitudes  of  each  trace  appear  the  same  size.  The  first  column  in  each 
figure  shows  the  analytic  solutions.  The  second  column  shows  the  same  seismograms 
band  pass  filtered  between  one  and  one  hundred  seconds.  The  third  column  shows  the 
bandpass  filtered  hybrid  solutions.  Figure  14  shows  the  results  using  the  0.15s  time 
function,  and  Figure  15  shows  the  results  for  the  0.45s  time  function.  Figure  16 
addresses  the  amplitude  correspondence  for  both  cases.  In  Figure  14  it  is  possible  to 
see  both  the  direct  arrival  and  the  arrival  reflected  from  the  free  surface  in  the  direct 
seismograms  recorded  at  depth.  These  arrivals  are  clearly  higher  frequency  than  the 
minimum  period  of  one  second  transmittable  by  the  FE  grid.  The  hybrid  solutions  for 
this  case  contain  unwanted  oscillations.  However,  when  the  hybrid  solutions,  and  the 
direct  solutions  are  both  band  pass  filtered  their  waveforms  are  in  excellent  agreement. 
Alternatively,  band  pass  filtering  the  forcing  functions  before  they  are  input  to  the  FE 
calculation  will  produce  unfiltered  hybrid  solutions  almost  indistinguishable  from  the 
filtered  hybrid  solutions  calculated  using  the  unfiltered  forcing  functions.  Aliasing  of 
energy  traveling  in  waves  of  frequency  higher  than  can  be  propagated  by  the  grid  into 
lower  frequency  bands  does  not  seem  to  be  a  problem  for  these  tests.  Figure  15  shows 
that  waveforms  of  the  direct  and  hybrid  solutions  correspond  well  when  the  input  has 
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Figure  15:  Effects  of  filtering  and  internal  damping  on  FE  solutions.  As  for  figure 
14  but  using  a  triangular  time  function  of  45s  duration. 


no  significant  high  frequency  component.  In  this  case  the  bandpass  filtering  makes  lit¬ 
tle  or  no  additional  improvement  to  the  correspondence  of  the  two  sets  of  results. 

The  spurious  oscillations  seen  in  the  unfiltered  hybrid  solutions  can  be  removed 
by  filtering  as  discussed  above,  or  by  the  introduction  of  a  damping  coefficient  within 
the  FE  calculations.  For  each  of  the  sets  of  forcing  functions  used  in  the  tests  dis¬ 
cussed  in  the  previous  paragraph  FE  solutions  were  determined  both  with  and  without 
damping.  The  damping  coefficient  used  in  the  FE  code  is  dependent  on  distance. 
Therefore,  a  single  damping  coefficient  cannot  produce  correct  amplitudes  at  a  variety 
of  distances.  At  the  distance  for  which  the  coefficient  is  calculated,  and  within  a  few 
grid  spacings  of  that  distance,  the  coefficient  produces  FE  results  that  correspond  very 
well  with  the  direct  solution.  If  the  direct  solution  contains  high  frequencies  the 
results  correspond  to  the  filtered  direct  results.  This  can  be  seen  by  examining  Figures 
14  and  15.  In  each  of  these  figures  column  two  shows  the  filtered  direct  solutions,  and 
column  four  shows  the  hybrid  solution  using  internal  damping  with  the  FE  calcula¬ 
tion.  The  damping  coefficient  is  chosen  to  be  optimal  at  the  illustrated  distance.  The 
correspondence  in  waveform  between  these  two  columns  is  very  good. 

Figure  16  shows  the  amplitude  ratios  for  the  tests  of  filtering  and  damping 
efficiency.  It  illustrates  the  relative  accuracy  of  using  the  damped  FE  and  the  filtered 
undamped  FE  approaches.  The  damped  FE  calculation  shows  more  variation  with 
depth  of  the  amplitude  ratios.  The  ratios  are  calculated  at  the  distance  for  which  the 
damping  coefficient  was  optimal.  However,  for  the  calculation  which  includes  damp¬ 
ing,  as  distance  from  this  point  increases  the  amplitude  ratios  rapidly  depart  from  one 
and  show  even  larger  variations  with  depth.  On  the  other  hand  the  ratios  for  the 
filtered  undamped  FE  calculation  remain  quite  stable. 

The  results  of  these  tests  can  be  summarized  by  noting  the  following  important 
observations.  If  the  input  seismograms  contain  frequencies  higher  than  those  that  can 
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Figure  18:  Amplitude  ratios  as  *  function  of  depth  for  seismograms  recorded  on 
a  depth  section  57  km  from  the  source. 


be  propagated  through  the  FE  grid,  as  in  the  case  of  the  rapid  rise  observed  for  the 
forcing  functions  and  analytic  seismograms  with  the  0.15s  time  function,  then  the  grid 
acts  as  a  filter.  The  amplitude  of  the  FE  solutions  are  significantly  reduced  with 
respect  to  the  analytic  solutions  and  spurious  oscillations  are  introduced  near  the  end 
of  the  waveform.  The  correspondence  of  amplitudes  and  waveforms  between  the  ana¬ 
lytic  solution  and  the  hybrid  solution  for  this  situation  is  poor  because  the  analytic 
solution  contains  higher  frequencies  than  the  hybrid  solution,  and  the  hybrid  solution 
contains  numerical  artifacts.  Low  pass  filtering  the  analytic  solution  and  the  hybrid 
solution  will  improve  the  waveform  fit  enormously,  and  bring  the  amplitude  ratios 
much  closer  to  one.  The  cutoff  frequency  of  the  low  pass  filter  should  be  at  most  the 
highest  frequency  that  the  grid  can  propagate.  It  is  not  important  whether  the  input 
forcing  functions  or  the  output  hybrid  solutions  are  filtered.  Both  approaches  produce 
equivalent  results.  When  input  forcing  functions  contain  no  frequencies  too  high  for 
propagation,  filtering  of  hybrid  solutions  or  analytic  synthetics  is  not  necessary  as  it 
does  not  change  the  goodness  of  fit.  The  use  of  internal  damping  in  the  FE  calcula¬ 
tion  is  distance  dependent  and  introduces  addition  variation  in  the  goodness  of  fit  of 
amplitude  and  waveform  as  a  function  of  depth.  Considering  all  these  observations 
the  best  combination  of  filtering  and  damping  when  a  range  of  distances  is  to  be  used 
is  to  filter  an  undamped  FE  solution.  If  theoretical  results  are  available  the  filtered 
undamped  solutions  should  be  compared  directly  to  them.  If  the  theoretical  results 
contain  high  frequencies  they  should  be  low  pass  filtered  before  being  compared  to  the 
filtered  undamped  solutions.  This  approach  should  give  the  best  correspondence  to 
the  theoretical  results. 

Teats  of  Accuracy:  SH  Lg  Mode  Sum  Input 

The  accuracy  of  coupling  a  2-D  SH  pulse  from  a  source  outside  a  FE  grid  into 
that  grid  have  now  been  thoroughly  discussed.  The  limitations  of  such  an  approach 
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have  also  been  defined.  However,  generalizing  these  results,  which  have  been  demon¬ 
strated  only  for  a  line  source  in  a  half-space,  to  a  line  source  or  a  point  source  in  a  lay¬ 
ered  media  or  even  in  a  layer  over  a  half-space  can  be  misleading.  Some  of  the  tests 
discussed  below  illustrate  cases  for  which  results  do  not  generalize.  The  physical 
model  used  in  all  of  these  tests  is  a  layer  over  a  half-space.  The  addition  of  the  layer 
above  the  half-space  allows  for  the  existence  of  Love  waves.  This  changes  the  behavior 
of  the  accuracy  of  coupling  a  distant  source  to  a  FE  calculation.  In  these  tests  the 
layer  has  a  thickness  of  thirty  two  kilometers,  an  S  wave  velocity  of  3.5  km/s,  and  a 
density  of  2.7  g/cc.  The  underlying  half-space  has  an  S  wave  velocity  of  4.5  km/s, 
and  a  density  of  3.5  g/cc.  The  initial  calculation  demonstrating  the  validity  of  cou¬ 
pling  a  distant  source  into  a  FE  grid  for  the  case  of  a  layer  over  a  half-space  used  a 
grid  spacing  of  0.53  km  in  both  x  and  z.  All  subsequent  tests  used  a  grid  spacing  of  0.5 
km  in  both  x  and  z.  All  tests  use  a  time  step  of  0.05  s.  Each  test  uses  the  same  set  of 
ninety  forcing  functions,  or  a  subset  of  that  same  set,  applied  to  the  leftmost  column 
of  nodes.  Forcing  functions  are  calculated  for  a  strike  slip  point  double  couple  source 
at  a  depth  of  eight  kilometers.  The  source  to  grid  edge  horizontal  distance  is  1500 
km. 

The  seismograms  used  as  forcing  functions  in  the  tests  discussed  below  are  calcu¬ 
lated  as  mode  sums  over  modes  with  periods  between  0.5  s  and  100  s.  The  fundamen¬ 
tal  mode  and  the  first  five  higher  modes  are  included  in  each  calculation.  This  is  not  a 
complete  representation,  but  produces  a  seismogram  that  shows  the  overall  character 
of  the  arrivals  seen  in  data  in  the  L,  group  velocity  range  of  between  2.8  and  3.5 
km/s.  A  complete  representation  would  require  as  much  as  five  times  the  amount  of 
calculation  to  produce  the  forcing  functions  or  the  synthetics  used  for  comparison  to 
the  hybrid  results.  To  assure  the  inclusion  of  all  energy  down  to  a  period  of  0.5  s 
twenty  two  higher  modes  would  be  necessary,  for  ail  energy  down  to  a  period  of  one 
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second  eleven  higher  modes  need  to  be  included.  The  cutoff  frequency  of  the  sixth 
higher  mode  is  1.9  s.  The  forcing  functions  used  are  complete  to  that  period  and  con¬ 
tain  a  significant  portion  of  the  energy  at  the  shorter  periods  considered.  For 
waveform  modeling  of  data  the  use  of  a  more  complete  mode  set  would  be  advisable, 
but  for  the  tests  discussed  here  it  is  an  unnecessary  expenditure  of  effort.  An  example 
of  a  forcing  function  or  direct  seismogram,  recorded  at  the  surface  at  a  distance  of 
1500  km  from  the  source  is  shown  in  Figure  17.  The  group  velocities  of  the  arrivals 
are  shown  along  the  horizontal  axis.  The  heavy  rectangle  enclosing  the  largest  portion 
of  the  seismogram  delimits  the  portion  of  the  trace  used  as  a  displacement  time  his¬ 
tory  to  drive  the  surface  node  on  the  depth  section  at  1500  km  from  the  source.  The 
seismogram  shown  in  this  figure  is  not  the  actual  trace  used  as  a  forcing  function.  It 
has  been  band  pass  filtered  between  0.01  and  1.0  Hz  and  has  had  the  WWSSN  short 
period  receiver  function  applied  to  it.  The  traces  used  as  forcing  functions  have  been 
band  pass  filtered  but  do  not  have  the  instrument  applied.  Examples  of  such  traces  at 
different  source  receiver  separations  are  shown  in  Figure  18.  In  this  figure  these  traces 
are  labeled  synthetic. 

At  this  point  it  is  useful  to  digress  and  to  more  precisely  define  the  RMS  ampli¬ 
tude  measure  used  in  subsequent  discussions  and  the  method  used  to  calculate  it.  The 
form  of  the  relation  used  to  determine  the  RMS  amplitude  is 


where  m  is  the  number  of  points  in  a  sampling  window  or  duration  Tr^s  seconds. 
The  location  of  the  sampling  window  with  respect  to  the  arrivals  of  maximum  ampli¬ 
tude  and  the  duration  of  the  trace  contained  within  the  sampling  window  have  a 
significant  effect  on  the  value  of  the  RMS  amplitude.  Care  must  be  taken  to  choose 
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Figure  17:  Sample  forcing  function  for  SH  L,  ware  FE  calculation.  Tbe  seismo¬ 
gram  includes  modes  with  periods  between  0.5s  and  1000s.  Only  the  fundamental  and 
the  first  five  higher  modes  are  used.  The  mode  sum  synthetic  has  been  transmitted 
1500  km  from  the  source  using  the  propagator  matrix  technique.  The  resulting 
seismogram  has  been  band  pass  filtered  between  0.01  and  1  hz,  and  has  the  WWSSN 
short  period  instrument  has  been  convolved  with  it.  The  portion  of  the  seismogram 
within  the  box  labeled  input  window  is  used  as  a  forcing  function.  The  numbers  along 
the  axis  indicate  the  group  velocities  of  the  arrivals. 


windows  for  two  sets  of  results  that  are  compatible  with  each  other  and  produce 
meaningful  comparisons.  For  example,  choosing  the  beginning  of  the  seismogram 
being  considered  as  the  first  point  in  the  sampling  window  can  cause  difficulties  when 
comparing  different  calculations  or  when  comparing  successive  seismograms  along  a 
depth  or  distance  section.  When  different  calculations  are  being  compared  extreme 
care  must  be  taken  to  insure  that  both  seismograms  begin  at  the  same  absolute  time. 
When  successive  seismograms  in  a  depth  or  distance  section  are  being  compared,  the 
duration  of  the  portion  of  the  trace  with  negligible  amplitude  that  occurs  before  the 
first  arrival  increases  with  distance  from  the  source.  This  creates  an  artifact  that 
makes  RMS  amplitudes  appear  to  decrease  with  distance.  To  avoid  these  and  other 
problems  some  other  method  of  determining  the  beginning  of  the  sampling  window 
must  be  devised.  It  would  be  possible  to  correct  for  the  distance  term  by  ignoring 
enough  points  at  the  beginning  of  the  seismogram  to  compensate  for  the  difference  in 
arrival  time.  However,  when  models  became  more  complicated  the  choice  of  the  first 
arrival,  and  the  determination  of  its  travel  time  from  the  source  becomes  more  compli¬ 
cated  due  to  the  possible  presence  of  strong  refracted  or  diffracted  phases.  Thus,  a 
more  general  approach  which  does  not  require  the  use  of  theoretical  travel  times  has 
been  chosen.  Two  different  methods  of  selecting  the  first  point  in  the  sampling  win¬ 
dow  have  been  used.  The  selection  of  one  of  these  two  methods  was  based  on  the 
characteristics  of  the  waveform  as  discussed  below. 

The  first  method,  which  was  used  in  the  majority  of  cases,  selects  the  first  point 
to  be  included  in  the  sampling  window  and  the  last  point  to  be  included  in  the  seismo¬ 
gram  by  bracketing  the  portion  of  the  seismogram  whose  RMS  amplitude  is  to  be 
measured.  This  piece  of  the  seismogram  is  bracketed  with  sections  of  seismogram  with 
sustained  low  amplitude  persisting  for  at  least  a  defined  minimum  duration.  This 
method  is  used  when  the  oscillations  in  the  waveform  are  very  small  until  the  first 


large  and  impulsive  arrival  occurs,  when  no  periods  of  sustained  low  amplitude  occur 
during  the  duration  of  the  trace,  and  when  the  amplitude  of  the  seismogram  decays  to 
a  sustained  low  amplitude  at  its  end.  Given  that  all  three  of  these  conditions  are 
satisfied  the  following  algorithm  produces  stable  and  accurate  RMS  amplitude  meas¬ 
ures.  The  algorithm  consists  of  five  steps: 

(1)  Determine  the  size  of  the  maximum  absolute  value  of  amplitude  in  the  seismo¬ 
gram 

(2)  Set  a  cutoff  value  for  the  amplitude  at  some  fraction  of  the  maximum  (usually 
.01  or  .05) 

(3)  Beginning  at  the  location  of  the  maximum  absolute  value  scan  toward  the  begin¬ 
ning  of  the  seismogram  until  a  series  of  samples  two  seconds  in  duration  all  have 
amplitudes  smaller  than  the  cutoff  amplitude.  The  first  point  in  this  series  is 
defined  to  be  the  first  point  in  the  sampling  window  used  to  determine  the  RMS 
amplitude. 

(4)  Beginning  at  the  location  of  the  maximum  absolute  value  scan  toward  the  end  of 
the  seismogram  until  a  series  of  samples  two  seconds  in  duration  all  have  ampli¬ 
tudes  less  than  the  cutoff  amplitude.  The  first  sample  in  that  series  of  samples 
defines  the  last  point  in  the  seismogram. 

(5)  Calculate  the  RMS  amplitude  for  sampling  windows  with  durations  15,  20,  25,  .., 
85s.  If  the  duration  of  the  sampling  window  exceeds  the  duration  of  the  seismo¬ 
gram  then  the  longest  sampling  window  used  is  the  duration  of  the  seismogram, 
and  that  duration  is  recorded. 

The  second  method  is  used  when  small  arrivals  precede  the  main  arrival  on  only 
a  subset  of  the  traces  being  considered,  when  the  small  arrivals  preceding  the  main 
arrival  contain  sustained  periods  of  low  amplitude  on  some  subset  of  the  traces,  and 
when  sustained  periods  of  low  amplitude  occur  elsewhere  within  any  of  the 


seismograms.  This  method  is  not  used  when  the  arrival  with  maximum  absolute 
amplitude  is  preceded  by  a  lower  amplitude  arrival  of  long  (>10s)  duration,  or  when 
the  maximum  amplitude  occurs  within  the  first  10s  of  the  seismogram.  This  method 
is  also  useful  when  slight  wraparound  occurs  as  a  result  of  convolving  the  instrument 
with  the  seismogram  producing  small  spurious  arrivals  before  the  main  peak.  If 
significant  wraparound  occurs,  the  convolution  length  is  extended  to  remove  it.  The 
determination  of  RMS  amplitude  is  not  expected  to  correct  it.  The  basics  for  the  algo* 
rithm  used  for  the  second  method  are  identical  to  those  used  in  the  first  method.  The 
only  differences  lie  in  the  method  used  to  select  the  first  point  in  the  sampling  window, 
and  in  the  duration  (5s)  of  sustained  low  amplitude  that  defines  the  end  of  the  seismo¬ 
gram.  When  using  the  second  method  the  first  point  in  the  sampling  window  is 
defined  to  be  the  point  in  the  time  series  for  a  time  ten  seconds  before  the  time  of  the 
peak  of  maximum  absolute  amplitude.  Clearly,  if  the  maximum  amplitude  occurs  late 
in  the  seismogram  this  method  is  not  useful.  However,  when  periods  of  sustained  low 
amplitude  occur  early  in  the  seismogram,  or  when  small  amplitude  early  arrivals  are 
detected  on  only  a  subset  of  the  seismograms  then  this  method  gives  stable  RMS 
amplitudes  whereas  the  amplitudes  using  the  first  method  will  not  be  stable. 

These  algorithms  yield  a  series  of  RMS  amplitude  values  that  can  be  used  to 
compare  results  derived  using  different  methods.  If  RMS  amplitudes  agree  within  all 
period  ranges,  and  waveforms  are  similar  the  fit  is  considered  to  be  excellent.  If 
discrepancies  occur  or  increase  with  the  length  of  the  sample  window,  then  there  is 
probably  a  problem  modeling  the  later  arrivals,  that  is,  the  representation  of  later 
arrivals  is  probably  incomplete.  If  the  shortest  time  windows  show  discrepancies  that 
decrease  rapidly  as  the  window  length  increases  then  the  amplitudes  of  the  initial 
arrivals  are  unstable  or  incorrect.  Examining  the  behavior  of  the  RMS  amplitude  as  a 
function  of  sample  window  length  can  give  insight  into  the  nature  of  and  the 
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underlying  reasons  for  misfit  between  methods.  In  most  of  the  discussions  that  follow, 
the  method  used  to  obtain  RMS  amplitudes  will  not  be  specified.  Both  methods 
described  above  give  consistent  interpretations  if  they  are  applied  to  all  seismograms 
in  a  given  example.  While  the  actual  amplitude  measures  may  vary  between  the  two 
methods  the  changes  seen  in  amplitude  within  a  given  example  are  the  same  regardless 
of  the  method  used  so  long  as  the  same  method  is  consistently  applied.  The  variation 
of  RMS  amplitude  with  the  length  of  the  sample  window  will  be  discussed  only  when 
such  variation  have  direct  impact  upon  the  results  of  the  comparisons. 

Now  I  will  return  to  the  discussion  of  the  RT  coupling  of  a  distant  line  source  in 
a  model  consisting  of  a  layer  over  a  half-space  into  a  FE  calculation.  First,  it  is  useful 
to  establish  that  the  coupling  of  energy  from  a  distant  source  into  a  FE  grid  is  valid 
for  a  layered  structure  more  complicated  than  a  half-space.  The  seismograms  illus¬ 
trated  in  Figures  18  and  19  verify  that  the  method  of  coupling  does  generalize  as  the 
theory  implies  it  should.  The  seismograms  in  Figure  18  are  recorded  at  two  distances 
from  the  grid  edge,  and  at  five  depths  for  each  of  these  distances.  A  pair  of  seismo¬ 
grams  is  shown  for  each  of  these  ten  locations.  The  upper  seismogram  in  each  pair  is 
the  direct  synthetic  determined  as  a  mode  sum  and  transmitted  to  the  receiver  using 
the  propagator  technique.  The  lower  seismogram  in  each  pair  is  the  hybrid  synthetic 
for  which  the  first  portion  of  the  path  is  traversed  using  the  propagator  matrix  tech¬ 
nique  and  the  remainder  of  the  path  is  traversed  in  a  FE  grid.  For  this  test  a  100x100 
node  grid  was  used  and  sixty  forcing  functions  were  applied  to  the  FE  grid  edge.  All 
seismograms  in  Figure  18  have  a  duration  of  fifty  five  seconds  and  have  been  band 
pass  filtered  between  .01  and  1  Hz.  The  amplitudes  of  all  traces  are  normalized  so  the 
maximum  peak  to  peak  amplitude  of  each  trace  appears  to  have  the  same  value. 
Within  each  pair  the  variation  in  peak  to  peak  amplitude  is  less  than  five  percent. 
The  peak  to  peak  amplitude  ratios  of  the  hybrid  synthetic  to  the  analytic  synthetic 
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Figure  18:  Comparison  of  hybrid  solution  to  mode  sum  synthetics.  The  upper 
truce  in  each  pair  is  u  mode  sum  synthetic  transmitted  using  only  propagator 
matrices.  The  lower  seismogram  in  each  pair  is  the  hybrid  solution.  The  depths  for 
each  row  and  the  horizontal  propagation  distances  in  the  FE  grid  for  each  column  are 
shown.  The  ratio  of  the  FE  to  synthetic  peak  to  peak  amplitudes  is  given  as  the 
number  below  each  pair  of  traces.  Each  trace  has  been  band  pass  filtered  between  .01 
and  1.0  Hz,  and  its  maximum  peak  to  peak  excursion  has  been  normalized  to  one. 


Figure  19:  Comparison  of  hybrid  solutions  and  3-D  mode  sum  synthetics.  The 
forcing  functions  used  are  3-D  mode  sum  synthetics.  The  upper  trace  in  each  pair  is 
the  mode  sum  synthetic,  the  lower  trace  is  the  hybrid  solution.  Numbers  at  the  left 
above  each  pair  of  seismograms  indicates  distance  from  source  to  receiver  at  the  sur¬ 
face.  Numbers  at  the  right  give  the  ratio  of  RMS  amplitude  of  the  hybrid  solution  to 
the  RMS  amplitude  of  the  direct  synthetic. 


are  shown  below  the  right  end  of  each  pair  of  traces.  The  major  differences  in  the 
synthetic  and  hybrid  solutions  are  seen  in  the  amplitude  of  the  first  large  positive  and 
negative  peaks.  The  amplitude  of  these  peaks  is  not  stable  since  they  are  the  result  of 
sampling  a  very  high  frequency  arrival  not  resolved  at  the  time  spacing  used.  These 
peaks  yield  the  maximum  peak  to  peak  amplitude  so  the  peak  to  peak  amplitude  is 
not  stable.  This  makes  the  peak  to  peak  amplitude  a  poor  measure  of  the  amplitude 
fit  between  models.  A  far  more  stable  measure  of  amplitude  correspondence  between 
seismograms  of  this  type  is  the  RMS  amplitude  calculated  over  some  time  window 
appropriate  to  the  seismograms  being  compared.  The  RMS  amplitudes  calculated 
with  window  lengths  of  as  much  as  45s  agree  to  better  than  one  percent  for  all  pairs  of 
seismograms  shown  in  Figure  18.  RMS  amplitudes  calculated  with  window  lengths 
greater  than  45  seconds  show  discrepancies  in  the  amplitude  correspondence  between 
each  pair  due  to  the  finite  length  of  the  seismogram.  Some  discrepancies  are  seen 
between  the  detailed  waveforms  of  the  direct  and  hybrid  synthetics.  These  discrepan- 
cies  are  of  at  least  two  types.  First,  and  most  easily  explained,  are  the  small  arrivals 
seen  at  the  beginning  of  the  seismogram  in  the  direct  synthetics  only.  The  lack  of 
these  arrivals  on  the  hybrid  synthetics  is  due  to  the  exclusion  of  the  very  low  ampli¬ 
tude  initial  portion  of  each  seismogram  used  as  forcing  function.  These  portions  lie 
outside  the  input  window  as  illustrated  in  Figure  17.  Second,  are  the  small  discrepan¬ 
cies  seen  in  the  higher  frequency  component  of  the  seismograms,  particularly  within 
the  first  twenty  seconds  of  the  trace.  To  better  quantify  these  shorter  period 
differences,  and  to  enable  one  to  view  the  record  as  one  would  see  data  from  an  event, 
the  WWSSN  short  period  instrument  response  is  convolved  with  seismograms  of  the 
type  shown  in  Figure  18  to  give  seismograms  of  the  type  shown  in  Figure  19. 

Figure  19  shows  comparisons  between  direct  and  hybrid  synthetics  for  a  selection 
of  horizontal  distances  along  the  free  surface.  Above  each  pair  of  seismograms  the 


horizontal  distance  from  the  source  to  the  receiver  is  shown.  For  this  test  a  FE  grid  of 
dimensions  150x300  nodes  was  used.  Beside  each  pair  of  traces  the  ratio  of  hybrid  to 
direct  RMS  amplitudes  is  given.  In  each  pair,  the  upper  seismogram  is  the  direct  syn¬ 
thetic  and  the  lower  trace  is  the  hybrid  synthetic.  The  duration  of  all  seismograms  is 
102s  and  the  peak  to  peak  amplitude  of  each  trace  is  scaled  to  the  same  value.  Each 
seismogram  shown  is  bandpass  filtered  between  0.01  and  1.00  Hz  and  includes  the 
WWSSN  instrument.  The  bandpass  filter  is  applied  to  the  forcing  functions  and  to 
the  resulting  FE  seismograms.  It  is  applied  twice  in  succession  to  the  direct  synthet¬ 
ics.  The  waveform  correspondence  between  the  two  types  of  synthetics  is  excellent  for 
short  horizontal  distances,  but  worsens  as  horizontal  distance  increases.  The  most 
striking  changes  occur  in  the  maximum  peak  to  peak  amplitudes,  and  in  the  ampli¬ 
tudes  of  successive  peaks  relative  to  each  other  in  the  first  ten  to  twenty  seconds  of 
the  record.  It  appears  that  the  waveform  correspondence  between  the  two  types  of 
synthetics  is  excellent  if  the  first  ten  to  twenty  seconds  of  the  record  is  omitted.  The 
discrepancies  early  in  the  seismogram  could  have  several  sources  such  as,  truncation 
errors  due  to  the  use  of  too  few  input  forcing  functions,  the  comparison  of  hybrid 
seismograms  which  experience  2-D  spreading  in  the  portion  of  their  paths  in  the  FE 
grid  to  direct  seismograms  which  experience  3-D  propagation  along  the  entire  path, 
reflections  from  the  bottom  of  the  FE  gTid,  and  possible  intrinsic  errors  due  to  the 
accuracy  of  the  numerical  codes  or  completeness  of  the  chosen  representation  for  the 
forcing  functions  and  direct  seismograms.  Some  of  these  possible  sources  for 
differences  in  waveform  will  be  discussed  in  detail  below  or  in  later  chapters.  It  is 
important  to  note  that  despite  the  small  variations  in  the  waveform  the  RMS  ampli¬ 
tude  correspondence  is  excellent  and  stable  as  propagation  distance  increases. 

Next  the  necessary  vertical  extent  of  input  forcing  functions  needed  to  produce 
an  accurate  mode  sum  hybrid  seismogram  and  the  truncation  error  associated  with 


that  vertical  extent  will  be  discussed.  For  the  case  of  the  SH  pulse  from  a  line  source 
discussed  earlier  it  was  found  that  accurate  results  required  a  vertical  extent  of  forcing 
functions  greater  than  the  depth  of  the  receiver  by  an  amount  that  increased  as  the 
horizontal  propagation  distance  increased.  This  meant  that  the  coupling  of  an  SH 
pulse  from  a  distant  source  into  a  grid  with  a  high  aspect  ratio  necessarily  produced 
results  that  became  very  poor  at  large  distances.  Fortunately,  this  behavior  does  not 
directly  generalize  to  the  case  of  a  layer  over  a  half-space.  For  a  layer  over  a  half¬ 
space  the  results  are  compatible  with  and  can  be  explained  by  considering  the  Love 
waves  to  be  modeled  in  terms  of  constructive  interference  of  post  critical  multiple 
reflections  trapped  within  the  crustal  layer.  This  interpretation  of  Love  waves  sug¬ 
gests  that  the  most  critical  nodes  at  which  forcing  functions  need  to  be  applied  are  the 
initial  nodes  in  each  row  of  the  grid  located  in  the  crustal  layer.  The  results  of  the 
calculations  illustrated  in  Figures  20  and  21  show  large  changes  in  the  waveforms  and 
amplitudes  of  hybrid  synthetics  occur  as  the  number  of  input  forcing  functions 
increases,  until  all  nodes  in  the  first  column  of  the  grid  that  lie  within  the  crustal  layer 
have  forcing  functions  applied  to  them.  Further  increases  in  the  depth  extent  of 
applied  forcing  functions  produce  only  small  changes  in  waveform  in  the  first  thirty 
seconds  of  the  trace  and  almost  no  changes  in  RMS  amplitudes.  This  remains  true  of 
the  hybrid  synthetics  after  propagation  distances  of  more  than  five  to  ten  layer 
thicknesses  of  the  crustal  layer. 

All  results  shown  in  Figures  20  and  21  were  calculated  in  a  FE  grid  of  dimension 
150x300  nodes.  Five  FE  calculations  were  completed,  each  using  a  different  number  of 
applied  forcing  functions.  These  tests  use  fifteen,  thirty,  forty  five,  sixty,  and  ninety 
applied  forcing  functions.  In  Figures  20  and  21  a  number  or  numbers  appear  below 
and  at  the  right  end  of  each  seismogram.  This  number  n  indicates  that  n  forcing 
functions  are  applied  to  the  topmost  n  nodes  in  the  first  column  of  the  FE  grid  in  the 


calculation  that  produced  the  illustrated  seismogram.  The  pairs  of  numbers  at  the  left 
end  of  each  seismogram  are  peak  to  peak  and  fifty  five  second  RMS  amplitude  ratios. 
For  each  seismogram  the  ratio  shown  is  the  ratio  of  the  amplitude  of  the  seismogram 
being  considered  to  the  amplitude  of  the  seismogram  at  the  same  location  for  the 
example  using  sixty  forcing  functions.  The  uppermost  number  is  the  peak  to  peak 
ratio,  the  number  below  it  is  the  RMS  ratio.  The  layer  thickness  of  thirty  two  kilome¬ 
ters  is  equivalent  to  sixty  four  rows  of  nodes.  The  sixty  fifth  row  of  nodes  is  the  first 
row  describing  the  half-space.  The  calculation  using  sixty  applied  forcing  functions  is 
used  as  a  reference  since  the  behavior  of  the  discrepancies  due  to  truncation  of  the 
vertical  extent  of  input  forcing  functions  changes  when  forcing  functions  are  applied 
outside  the  layer. 

Figure  20  shows  three  groups  of  seismograms.  Each  group  is  recorded  at  the  sur¬ 
face  at  a  given  horuontal  distance.  That  distance  is  shown  on  the  figure  above  each 
group  of  traces.  The  first  group  consists  of  a  single  seismogram  recorded  five  kilome¬ 
ters  or  ten  elements  from  the  grid  edge.  This  seismogram  represents  the  identical 
results  observed  in  the  five  test  cases.  The  second  group  of  seismograms  are  recorded 
twenty  five  kilometers  from  the  grid  edge.  The  first  seismogram  in  this  group 
represents  the  almost  identical  results  for  the  three  test  cases  using  the  largest 
numbers  of  applied  forcing  functions.  Amplitude  discrepancies  between  these  cases 
were  less  that  0.2%,  and  no  differences  in  waveform  could  be  seen  when  the  seismo¬ 
grams  for  each  case  were  overlaid.  However,  the  second  two  seismograms  in  the 
group,  representing  the  tests  with  thirty  and  fifteen  applied  forcing  functions  respec¬ 
tively  show  significant  differences  in  waveform  and  amplitude.  The  most  striking 
effect  of  using  fewer  applied  forcing  functions  is  the  reduction  in  the  amplitudes  of  the 
later  arrivals  in  the  seismograms  with  respect  to  the  early  arrivals.  The  peak  to  peak 
amplitudes  and  the  RMS  amplitudes  with  longer  window  lengths  were  also  reduced  as 


Figure  20:  Comparison  of  hybrid  results  using  different  numbers  of  input  forcing 
functions.  Each  group  of  seismograms  shows  results  at  a  given  horizontal  distance,  as 
labeled  above  that  group.  The  number  or  numbers  below  and  to  the  right  of  each 
seismogram  indicate  the  number  of  nodes  where  forcing  functions  are  applied.  The 
pairs  of  numbers  to  the  left  of  each  seismogram  are  amplitude  ratios.  The  upper 
number  in  each  pair  is  the  ratio  of  the  peak  to  peak  amplitudes  of  that  trace  to  the 
trace  using  60  forcing  functions  recorded  at  that  distance.  The  tower  number  in  each 
pair  is  the  corresponding  ratio  for  RMS  amplitudes. 


the  number  of  applied  forcing  functions  decreased.  In  particular,  for  the  seismogram 
generated  using  thirty  applied  forcing  functions  the  RMS  amplitude  measure 
corresponds  well  with  the  examples  using  more  forcing  functions  until  the  window 
length  of  the  RMS  operator  exceeds  sixty  seconds,  then  the  discrepancy  increases  as 
the  window  length  increases.  For  the  example  using  fifteen  forcing  functions  the  trend 
seen  in  the  RMS  amplitudes  follows  the  same  pattern  but  the  differences  in  peak  to 
peak  amplitude  mean  that  agreement  is  not  good  even  at  a  window  length  of  fifteen 
seconds.  If  results  of  the  calculation  using  fifteen  forcing  functions  are  examined  at 
locations  between  the  two  distances  illustrated  it  is  found  that  the  behavior  seen  in 
the  thirty  forcing  function  example  is  seen  at  distances  around  twelve  kilometers  in 
the  fifteen  forcing  function  calculation.  The  third  group  of  seismograms  is  recorded  at 
a  distance  of  fifty  kilometers.  The  second  two  seismograms  in  that  group  also  show 
the  rapid  decline  in  RMS  amplitude  as  the  RMS  window  length  increases.  However, 
the  example  using  thirty  forcing  functions  now  also  shows  poor  agreement  even  in  the 
fifteen  second  window.  The  example  using  forty  five  forcing  functions  is  not  illus¬ 
trated  at  this  distance.  The  RMS  amplitude  of  the  forty  five  forcing  function  seismo¬ 
gram  corresponds  well  with  the  more  complete  examples  for  RMS  windows  as  long  as 
sixty  five  seconds.  It  was  not  considered  necessary  to  illustrate  this  seismogram  as  it 
was  so  similar  to  the  first  seismogram  in  the  group.  From  examining  the  seismograms 
at  these  three  distances  patterns  are  beginning  to  emerge,  that  will  be  clarified  by  the 
additional  distances  illustrated  in  the  next  figure. 

Figure  21  is  a  continuation  of  Figure  20  for  seismograms  recorded  at  greater  dis¬ 
tances  along  the  free  surface.  The  amplitude  ratios  of  the  two  examples  using  fifteen 
and  thirty  forcing  functions  continue  to  decrease  with  distance,  although  the  RMS 
amplitudes  have  almost  stabiized  by  the  time  the  horizontal  distance  has  reach  150 
km.  The  disturbances  visible  in  the  trace  have  become  larger  with  respect  to  the 


shrinking  arrival  that  determines  the  maximum  peak  to  peak  amplitude  of  the  seismo¬ 
gram.  For  the  example  using  forty  five  forcing  functions  the  later  arrivals  continue  to 
shrink  with  respect  to  the  earlier  arrivals  as  horizontal  distance  increases,  reducing  the 
RMS  amplitude  in  successively  shorter  time  windows.  At  distances  between  seventy 
five  and  one  hundred  kilometers  the  RMS  amplitude  stabilizes  at  about  90%  of  the 
values  seen  for  the  examples  using  more  forcing  functions.  However,  the  waveforms 
correspondence  continues  to  decay.  Generalizing  the  results  seen  in  the  examples 
using  a  depth  extent  of  forcing  functions  considerably  less  than  the  thickness  of  the 
crustal  layer  leads  to  several  observations.  Seismograms  along  the  surface  are  identical 
regardless  of  the  number  of  forcing  functions  used  provided  that  the  propagation  dis¬ 
tance  in  the  FE  grid  is  less  than  or  equal  to  the  depth  extent  of  the  forcing  functions. 
Beyond  that  distance  range  the  amplitude  of  the  waveform  begins  to  decay.  The 
amplitudes  of  the  latest  arrivals  are  the  first  to  be  reduced.  As  the  distance  increases 
progressively  earlier  arrivals  are  reduced  in  amplitude  while  amplitudes  of  the  later 
arrivals  continue  to  decay  further.  As  the  number  of  forcing  functions  increases  the 
rate  of  amplitude  decay  decreases  and  the  RMS  amplitude  value  at  which  the  ampli¬ 
tude  eventually  stabilizes  increases.  When  the  depth  extent  of  forcing  functions 
exceeds  the  crustal  layer  thickness  this  behavior  changes.  At  a  distance  of  three  hun¬ 
dred  kilometers  the  two  examples  using  sixty  and  ninety  forcing  functions  continue  to 
agree  very  well.  They  show  some  differences  in  waveform,  barely  perceptible  when  the 
waveforms  are  overlaid,  in  the  relative  amplitudes  of  successive  peaks  during  the  first 
twenty  to  thirty  seconds  of  the  seismogram.  At  larger  distances  this  discrepancy 
increases  slightly,  but  the  amplitude  ratios  remain  very  stable  and  the  waveform 
differences  remain  small  and  confined  to  the  first  thirty  seconds  of  the  seismograms. 
Thus,  truncation  error  can  explain  at  least  a  part  of  the  waveform  discrepancies  seen 
in  the  early  portions  of  the  seismograms  illustrated  in  the  surface  section  shown  in 


Figures  22  and  23  illustrate  the  effects  of  using  different  numbers  of  row  of  nodes, 
that  is  varying  the  vertical  extent  of  the  FE  grid,  on  the  resulting  SH  L,  mode  sum 
seismograms.  In  each  of  these  figures  two  groups  of  three  seismograms  are  shown. 
Each  group  is  labeled  with  its  propagation  distance  in  the  FE  grid.  The  same  forcing 
functions,  for  a  source  1500  km  distant  at  a  depth  of  eight  kilometers,  are  used  for 
each  calculation.  All  seismograms  in  these  figures  are  plotted  at  the  same  scale.  The 
uppermost  seismogram  in  each  group  was  calculated  in  a  FE  grid  with  ninety  rows, 
the  center  seismogram  in  a  FE  grid  with  one  hundred  fifty  rows  and  the  lower  seismo¬ 
gram  in  a  FE  grid  with  three  hundred  rows.  An  additional  difference  between  the 
grids  used  in  these  tests  is  that  the  grid  in  which  the  uppermost  seismogram  in  each 
group  is  calculated  has  a  length  of  five  hundred  nodes  and  the  other  two  grids  have  a 
length  of  three  hundred  nodes.  Consequently,  differences  seen  in  the  last  ten  to  twelve 
seconds  of  the  upper  two  records  can  be  attributed  to  reflections  from  the  righthand 
edge  of  the  grid.  All  other  differences  seen  between  these  seismograms  can  be  attri¬ 
buted  to  reflections  from  the  bottom  edge  of  the  grid.  Since  the  travel  time  to  the 
bottom  of  the  grid  is  increased  as  the  number  of  rows  in  the  grid  is  increased,  any 
reflections  from  the  bottom  of  the  grid  should  be  seen  at  different  points  in  the  seismo¬ 
gram  for  each  grid.  The  length  of  the  seismogram  is  chosen  so  that  reflections  from 
the  bottom  of  the  three  hundred  row  grid  should  arrive  after  the  end  of  the  seismo¬ 
gram.  The  reflections  from  the  righthand  edge  of  the  grid  are  small  but  clearly  visible, 
and  are,  as  expected,  identical  on  both  of  the  lower  two  seismograms  in  each  group. 
Examining  the  group  of  seismograms  for  a  propagation  distance  of  fifty  kilometers 
within  the  FE  grid  shows  essentially  no  differences  between  the  seismograms  recorded 
in  grids  with  ninety  and  three  hundred  rows.  When  compared  carefully  to  these 
seismograms  the  seismogram  calculated  in  the  grid  with  one  hundred  fifty  rows  shows 
some  small  differences  in  the  amplitudes  of  some  of  the  peaks.  The  RMS  amplitudes 
agree  to  within  less  than  0.1%.  This  implies  that  reflections  from  the  bottom  of  the 


Figure  22:  Comparison  of  hybrid  solutions  for  cases  using  different  numbers  of 
rows  of  nodes.  Results  show  wide  angle  reflections  are  not  important.  Each  group  of 
three  seismograms  is  recorded  at  the  free  surface  at  the  indicated  distance  from  the 
grid  edge.  The  upper  trace  in  each  group  uses  90  rows  of  nodes,  the  center  trace  150 
rows  of  nodes  and  the  bottom  trace  300  rows.  The  peak  to  peak  amplitudes  of  each 
trace  are  normalized  to  one.  Differences  in  peak  to  peak  amplitudes  are  so  small  that 
all  peak  to  peak  and  RMS  amplitude  ratios  are  1.00.  Duration  of  each  seismogram  is 
102  seconds. 
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grid  are  insignificant  with  short  path  lengths  in  the  grid.  For  a  path  length  of  one 
hundred  kilometers  within  the  FE  grid  the  reflections  from  the  bottom  of  the  grid 
remain  insignificant,  as  illustrated  in  the  second  group  of  seismograms  in  figure  22. 
The  two  groups  of  seismograms  in  Figure  23  show  that  the  reflections  from  the  bot¬ 
tom  of  the  grid  remain  insignificant  to  distances  of  at  least  one  hundred  fifty  kilome¬ 
ters.  Thus,  even  for  aspect  ratios  as  high  as  three  the  reflections  from  the  bottom  of 
the  grid  are  not  significant  when  the  propagation  of  SH  type  L,  mode  sum  seismo¬ 
grams  in  a  layer  over  a  half-space  is  being  considered.  This  conclusion  will  be  further 
supported  when  the  transparent  boundary  conditions  are  considered  in  the  following 
chapter.  It  will  be  demonstrated  that  reflections  from  the  bottom  of  the  grid  are 
significant  only  when  structures  are  not  plane  layered. 

Summary 

In  this  chapter  the  basics  of  the  modal  propagator  matrix  method  and  the  FE 
method  were  discussed,  including  the  modifications  made  to  the  existing  implementa¬ 
tions  of  these  methods.  An  approach  to  apply  the  Representation  Theorem  to  com¬ 
bine  these  techniques  into  a  hybrid  method  was  explained.  The  particulars  of  the 
implementiation  of  the  coupling  technique  differ  depending  on  whether  the  propagator 
results  are  to  be  passed  into  the  FE  calculation  or  the  the  FE  results  are  to  be 
transmitted  further  using  the  propagator  matrix  technique.  The  coupling  from  modal 
propagator  matrix  results  to  the  FE  method  was  demonstrated  for  the  simple  example 
of  a  layer  over  a  half-space.  This  simple  structure  allowed  the  calculation  of  analytic 
synthetics  to  which  the  hybrid  synthetics  could  be  directly  compared  to  give  estimates 
of  the  accuracy  of  the  coupling  technique. 
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